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ABSTRACT 


This  research  compares  the  performance  of  three  filters 
which  have  been  applied  to  the  problem  of  orbit  determina¬ 
tion  using  actual  satellite  tracking  data  obtained  from 
ground  based  radars.  The  states  estimated  are  the  osculat¬ 
ing  classical  orbital  elements  and  the  satellite  ballistic 
coefficient.  The  dynamics  used  to  propagate  the  state 
vector  forward  include  the  two-body  acceleration  plus 
perturbations  due  to  atmospheric  drag,  zonal  harmonics  in 
the  geopotential  through  and  tesseral  harmonics  in  the 
geopotential  through  The  atmospheric  density  model 

used  is  an  exponential  model  that  includes  diurnal  varia¬ 
tions  and  variations  in  the  decimeter  solar  flux.  The 


observations  used  to  update  the  state  vector  estimates  are 
slant  range,  azimuth,  and  elevation  relative  to  a  radar 


site. 


The  three  filters  investigated  in  this  research  are  a 


nonlinear  least  squares  filter,  an  Extended  Kalman  filter, 
and  a  Gauss  second  order  filter.  Data  are  processed  for 


three  different  satellites.  The  first  is  a  high  altitude 


(1000km  at  perigee),  non-circular  {e=0.015),  orbit.  The 


second  satellite  orbit  is  a  low  altitude  (250km  at  perigee), 
non-circular  (e*0.01),  orbit.  The  final  orbit  is  a  low 
altitude  (300km),  nearly  circular  (e=0.0003),  orbit. 


The  filters  are  compared  using  four  criteria:  estima¬ 
tion  errors,  prediction  errors,  computer  time  of  operation, 
and  computer  storage  requirements.  The  Gauss  second  order 
filter  is  shown  to  provide  a  substantial  improvement  in 
orbit  determination  accuracy  for  satellites  subject  to 
significant  perturbing  accelerations. 
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CHAPTER  1 


INTRODUCTION 

1.1.  Problem  Statement 

The  orbit  determination  problem  consists  of  two  basic 
parts;  propagation  of  the  state  estimates  forward  in  time 
and  updating  the  state  estimate  based  upon  new  measurements 
of  parameters  which  are  functions  of  the  states.  The  prob¬ 
lem  is  complicated  by  the  following  three  factors.  First, 
the  initial  state  vector,  the  satellite's  orbital  elements, 
is  not  known  exactly.  Next,  the  dynamics  model  used  to 
propagate  the  state  forward  is  only  an  approximate  model. 
Finally,  the  measurements  are  corrupted  by  noise. 

The  states  to  be  estimated  are  the  classical  orbital 
elements  and  the  ballistic  coefficient  of  the  satellite. 
These  elements  are: 

a  -  semi-major  axis, 
e  -  eccentricity, 
i  -  inclination, 

n  -  longitude  of  ascending  node, 
oj  -  argument  of  perigee, 

M  -  mean  anomaly,  and 
B  -  ballistic  coefficient 

The  classical  orbital  elements  were  chosen  due  to  their 
descriptive  nature.  The  semi-major  axis  and  eccentricity 


determine  the  size  and  shape  of  the  orbit.  The  inclina¬ 
tion  and  longitude  of  ascending  node  describe  the  orien¬ 
tation  of  the  orbital  plane  in  inertial  space.  The 
argximent  of  perigee  describes  the  orientation  of  the  orbit 
in  the  orbital  plane.  Finally,  the  mean  anomaly  locates 
the  satellite  within  the  orbit.  The  ballistic  coefficient 
is  also  estimated  since  it  is  not  known  for  many  objects 
orbiting  the  Earth  and  its  estimate  can  incorporate  uncer¬ 
tainty  in  the  atmospheric  model. 

The  dynamics  model  used  to  propagate  the  state  vector 
forward  in  time  for  a  satellite  must  contain  the  two-body 
equations  of  motion  and  any  perturbations  which  are  con¬ 
sidered  significant.  The  method  used  in  this  work  is  a 
variation  of  parameters  formulation  of  the  equations  of 
motion  which  are  integrated  using  special  perturbations 
theory.  The  perturbations  include  atmospheric  drag,  zonal 
terms  J2  ~  Jg  in  the  geopotential,  and  tesseral  terms 
J22  “  J44  in  the  geopotential.  Since  no  model  is  exact, 
there  are  errors  in  the  equations  of  motion.  These  errors, 
called  process  noise,  are  modeled  as  zero-mean,  white, 
Gaussian  noise. 

Similarly,  the  measurements  used  to  update  the  states 
are  not  exact,  but  are  corrupted  by  noise  in  the  observation 
data.  These  errors  are  also  modeled  as  zero-mean,  white. 


Gaussian  noise. 


The  problem  of  satellite  orbit  determination  is  a 
very  important  one  which  has  received  much  attention  in 
the  literature.  In  this  work  we  will  compare  the  accuracy 
and  efficiency  of  several  different  filters  for  orbit 
determination-  Among  previous  researchers  who  have  made 
such  comparisons  are  Kolenkiewicz  and  Fuchs  (1)  and  Myers 
(2) .  The  standard  for  comparison  is  a  batch,  nonlinear, 
least  squares  filter.  This  method  has  been  used  for  years 
and  has  proven  to  be  accurate  when  applied  to  the  problem 
of  orbit  determination  for  a  large  number  of  satellites 
in  different  orbits.  Currently,  the  Air  Force  Space  Com¬ 
mand  uses  a  nonlinear  least  squares  filter  for  element  set 
maintenance  for  all  objects  orbiting  the  Earth  (3) .  These 
objects,  including  active  payloads,  rocket  bodies,  and 
satellite  fragments,  number  over  5000.  The  other  filters 
considered  are  all  recursive  filters.  Shavers  (4)  compares 
the  performance  of  batch  and  recursive  filters  for  orbit 
determination.  Among  the  advantages  of  recursive  filtering 
are  a  state  estimate  based  on  the  most  recent  data  avail¬ 
able  and  the  lack  of  a  need  to  store  the  data  for  batch 
processing.  One  type  of  recursive  filter  is  the  Extended 
Kalman  Filter.  The  final  type  of  filter  is  a  second  order 
filter,  which  derives  its  name  from  the  inclusion  of  the 
second  order  terms  from  the  Taylor  series  expansion  used 
to  linearize  the  filter  dynamics  and  observation  equations. 


The  second  order  terms  occur  in  the  dynamics,  observations, 
and  filter  gain  equations.  Tapley  and  Choe  (5)  compared 
the  performance  of  a  second  order  filter  with  the  perform¬ 
ance  of  an  Extended  Kalman  filter  for  interplanetary  orbit 
determination.  They  investigated  the  effect  of  the  second 
order  terms  in  the  dynamics,  observations,  and  filter 
gain  separately,  concluding  that  the  second  order  term  in 
the  filter  gain  resulted  in  the  most  significant  improvement 
over  the  Extended  Kalman  filter.  Taylor  (6)  calculated 
the  magnitude  of  the  bias  terms  in  the  dynamics  and  obser¬ 
vations  for  the  determination  of  an  orbit  of  an  Earth 
satellite.  He  concluded  that  inclusion  of  the  observation 
bias  terms  should  improve  the  filter  performance,  but  that 
the  inclusion  of  the  dynamic  bias  terms  should  not  alter 
filter  performance.  However,  Taylor  (6)  did  not  verify 
or  quantify  the  improvement  due  to  second  order  terms  as 
they  were  never  employed  in  the  filters  he  used.  Athens, 
Wishner,  and  Bertolini  (7)  compared  the  performance  of  a 
second  order  filter  to  the  perfoxnnance  of  a  first  order 
filter  when  applied  to  the  problem  of  estimating  the  posi¬ 
tion  and  velocity  of  a  vertically  falling  body.  They  con¬ 
cluded  that  the  second  order  filter  yields  superior  perform¬ 
ance  when  nonlinearities  are  significant.  They  also 
concluded  that  the  major  improvement  was  due  to  the  bias 
correction  term  in  the  dynamics. 


In  this  research  the  orbit  determination  accuracy  and 
efficiency  of  two  recursive  filters,  an  iterated-extended 


Kalman  filter  (EKF)  and  a  Gaussian  second  order  filter  (GSF) 

will  be  compared  to  that  of  a  batch  least  squares  filter 

(LSF) .  Actual  observation  data  for  satellites  in  three 

different  types  of  orbits  will  be  used.  The  Gaussian 

second  order  filter  will  be  examined  in  greatest  detail. 

* 

Following  the  approach  of  Tapley  and  Choe  (5),  the  effect 
of  the  second  order  terms  in  the  dynamics,  observations, 
and  gain  will  be  evaluated  individually.  The  efficiency  of 
the  filters,  expressed  as  computer  time  of  operation  and 
required  storage,  will  be  compared. 

1 . 2  Assumptions 

The  assumptions  made  in  this  research  are: 

a.  The  system  dynamics  are  continuous  with  proces 
noise  modeled  by  zero-mean,  white,  Gaussian  noise. 

b.  The  measurements  are  discrete  and  are  cor¬ 
rupted  by  zero-mean,  white,  Gaussian  noise. 

c.  The  process  noise  and  measurement  noise  are 
uncorrelated . 

d.  The  satellite  ballistic  coefficient  is  an 
unknown  constant. 

e.  The  most  significant  perturbing  forces  are 
due  to  atmospheric  drag  and  the  geopotential  terms 


J2  -  Jg  and  J22  "  '^44*  other  forces  may  be  neglected. 

f.  The  satellites  are  considered  to  be  non¬ 
maneuvering. 

1.3.  Thesis  Overview 

The  purpose  of  this  thesis  is  to  compare  the  perform¬ 
ance  of  various  filters  applied  to  the  orbit  determination 
problem.  Chapter  2  is  a  description  of  the  mathanatical 
system  dynamics  model  used  to  propagate  the  state  vector 
forward  and  the  observations  used  to  update  the  estimates 
of  the  states. 

Chapter  3  provides  an  overview  of  the  estimation  prob¬ 
lem  and  a  detailed  description  of  each  of  the  filter 
algorithms  used  in  this  research. 

Chapter  4  describes  the  satellite  orbits  and  the  data 
used  to  compare  the  filter  performance. 

Chapter  5  presents  the  results  obtained  by  the  filters 
using  actual  observational  data.  Results  are  presented 
regarding  estimation  errors,  prediction  errors,  computer 
storage  requirements,  and  computer  operation  time. 

Finally,  Chapter  6  contains  the  conclusions  of  this 
research  and  recommendations  for  future  work. 


CHAPTER  2 


SYSTEM  MODEL 


2.1.  Introduction 

The  propagation  of  the  satellite  state  vector  requires 
an  accurate  force  model  in  the  equations  of  motion  and  a 
means  of  integrating  the  equations  of  motion.  The  method 
used  in  this  research  is  to  express  the  equations  of  motion 
in  terms  of  Gauss'  variation  of  parameters  equations.  The 
equations  of  motion  are  then  integrated  nximerically  using 
the  method  of  special  perturbations. 

In  addition  to  the  dynamics  model,  the  mathematical 
relationship  between  the  observations  and  states  must  be 
defined  in  order  for  the  filters  to  process  the  data  and 
update  the  estimates  of  the  states.  The  observations  used 
in  this  research  are  slant  range  (p),  azimuth  (Az) ,  and 
elevation  (El)  relative  to  the  radar  site. 


2.2.  Dynamics 

The  equations  of  motion  of  a  satellite  orbiting  the 
Earth  are  given  by  Newton's  2^*^  Law  as: 


where: 


r  + 


2.1 


r  is  the  position  vector  of  the  satellite  from 


the  center  of  the  Earth, 


8 


U  is  the  Earth's  gravitational  constant,  and 
ap  are  the  perturbing  accelerations. 

These  equations  can  be  integrated  directly,  an  approach 
called  Cowell's  method,  but  this  approach  is  rather  slow 
and  inefficient  due  to  the  large  changes  in  the  magnitudes 
of  the  components  of  the  satellite  position  vector  and 
velocity  vector.  An  alternate  approach,  called  the  Varia¬ 
tion  of  Parameters,  takes  advantage  of  the  periodic  nature 
of  equation  2.1.  If  only  the  two-body  system  is  considered, 
equation  2.1  becomes: 

r  +  r  =  0.  2.2 

_  r^  -  - 

If  the  equations  of  motion  are  now  written  in  terms  of  the 
classical  orbital  elements,  a  system  of  six  first  order 
differential  equations  result  in  which  all  the  derivatives 
are  ecjual  to  zero  except  for  mean  anomaly  which  is; 

=  n,  2.3 


wher  e : 

n  =  /u/a^ ,  the  mean  motion. 

Writing  the  perturbing  accelerations  in  terms  of  the  classi¬ 
cal  orbital  elements,  the  system  equations  of  motion  become: 

k  =  f (x(t) ,t)  2.4 
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where ; 

X  is  the  seven-component  state  vector  containing 
the  classical  orbital  elements  and  the  satellite  ballistic 
coefficient,  and 

£  is  the  vector  containing  the  perturbing  accel¬ 
erations  and  the  mean  anomaly  rate  (Eq.  2.3) . 

The  perturbations  most  significant  for  low-Earth  orbit 
determination  are  acceleration  due  to  atmospheric  drag  and 
acceleration  due  to  a  non-spher ical  Earth.  The  mathematical 
model  of  acceleration  due  to  atmospheric  drag  is: 

2.5 

where : 

a^  is  the  acceleration  due  to  atmospheric  drag, 

B  is  the  satellite  ballistic  coefficient, 
p  is  the  atmospheric  density,  and. 
v_  is  the  velocity  of  the  satellite  relative  to 
the  atmosphere. 

The  accelerr ♦•ion  due  to  atmospheric  drag  must  now  be  ex¬ 
pressed  in  terms  of  the  state  vector,  x(t) .  The  approach 
used  in  this  thesis  is  that  developed  by  King-Hele  (8) . 

The  resulting  variational  equations  for  an  oblate  atmos¬ 
phere  with  rotational  velocity  equal  to  the  Earth's 
rotational  velocity  are: 
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a  =  -pBFa^v^/u 
k  =  pBFv(e+cos(f) ) 

i  =  -pvr^ojg  (F/up)  sin  ( i)  cos^  (u/2)  2.6 

2  1/2 

=  pvr  ojg(F/up)  '  sin  (u)  cos  (u) 
d)  =  -  cos(i) 

where; 

F  =  (1-r  0)  cos(i)/v  ) 

P  “  P 

Tp  is  the  radius  of  perigee, 

V  is  the  velocity  at  perigee, 

P 

f  is  the  satellite  true  anomaly, 
u  =  (ii+f, 

(Ug  is  the  Earth's  rotational  velocity,  and 
p  is  the  semi-latus  rectum. 

A  major  computational  effort  in  calculating  the  per¬ 
turbation  due  to  atmospheric  drag  is  in  determining  the 
density,  p.  The  simplest  atmospheric  model  assumes  a 
spherically  symmetric  density  which  varies  exponentially 
with  altitude; 

p  =  p^exp(-(h-h^) /H)  2.7 

wher  e ; 

p^  is  the  density  at  h^, 

h  is  the  altitude  above  the  surface  of  the  Earth, 

h^  is  the  reference  altitude,  and 
o 

H  is  the  scale  height. 


A  more  sophisticated  model  developed  by  Jacchia  (9) 
includes  a  diurnal  variation  due  to  solar  heating  and 
variations  in  the  decimeter  solar  flux,  F10.7.  The  model 
used  in  this  research  is  the  1960  Jacchia  model  with 
modifications  described  in  Reference  3.  For  this  model  the 
density  for  altitudes  less  than  700km  is; 

p  =  p^0.85F10.7[l+0.2375(exp(0.0102h)-1.9)  (1+cos  (ij; )  )  ^] 

2.8 

where; 

logioP^  =  -15.738-0. 00368h+6. 363exp(-0. 0048h) , 
p  is  the  density  in  slug/ft^, 
Fl0.7=1.5+0.8cos(iTd/2010)  , 
d  is  the  number  of  days  since  31  Dec.  1957, 
h  is  the  altitude  above  the  surface  of  the  Earth 
in  nautical  miles,  and 

ij;  is  the  angle  between  the  satellite  and  the 
subsolar  atmospheric  bulge. 

Above  700km  the  density  is 

p  =  0. 00504F10.7h"®  [1/8  (l+cos;Jj)  ^  (h^-6xl0^) +6x10®]  2.9 

The  density  scale  height,  H,  varies  with  altitude  and  is 
obtained  from  a  stored  table. 

The  other  major  perturbation  affecting  the  motion  of 
a  satellite  in  a  low-Earth  orbit  is  caused  by  the 


12 


gravitational  attraction  of  the  non-spherical  Earth.  The 
perturbing  potential  due  to  the  axially  symmetric  part 
of  the  Earth's  gravitational  field  is  (10): 

00 

F=(u/r)  Z  J„(R/r)^P  (sin(L))  2.10 

n=2  "" 

where: 

the  are  the  zonal  harmonic  coefficients  of  the 
Earth's  mass  distribution, 

R  is  the  Earth's  mean  equatorial  radius, 

are  Legendre  polynomials  of  degree  n,  and 
L  is  the  geocentric  latitude,  sinL=sinisin (oi+f )  . 
The  perturbing  acceleration  may  be  determined  from  the 
potential  function: 

ap  =  VF.  2.11 

A  convenient  coordinate  frame  for  expressing  the  compo¬ 
nents  of  acceleration  is  the  STW  frame  where: 

S  is  the  component  along  the  radius  vector, 

T  is  perpendicular  to  S  in  the  orbit  plane  in 
the  direction  of  motion,  and 

W  is  normal  to  the  orbit  frame  so  that  the  STW 
frame  forms  a  right-hand  coordinate  set. 

The  components  of  the  disturbing  acceleration  are  then: 

S  =  3F/3r, 

T  =  (l/r)3F/3u,  and  2.12 

W  =  (l/rsin(u) ) 3F/3i. 


The  components  of  the  disturbing  acceleration  due  to 


! 


} 


the  first  five  terms  (n=2,...,6)  of  the  disturbing  geopoten¬ 
tial,  eq.  2.9,  are  given  by  Merson  (11): 

S=(u/r^) [ (3/2) J2(R/r) ^(3A^-1)+2J3 (R/r) ^ (5A^-3A) 

+  (3/8)  (R/r)  ^  (3  5a'^-3  0A^+3) 

+(3/4) J5(R/r) ^(63A^-70A^+15A) 

+  (7/16)  Jg(R/r)  ^(231A®-315a'^+105A^-5)  ]  ,  2.13 

T=-(y/r^)  sinicosu[3J2(R/r)  ^A+(3/2)  J3  (R/r)  ^  (5A^-1) 
+(5/2) (R/r) ^(7A^-3A) 

+  (15/8)  J5  (R/r)  ^(21a'^-14A^+1) 

+(21/8) Jg (R/r)® (33A^-30A^+5A] ,  2.14 

2 

W=(u/r  ) cosi [bracketed  term  in  eq.  2.14).  2.15 

where : 

A  =  sin(i)sin(u) 

J2  =  1082.628x10”® 

J3  =  -2.538xl0"® 

J4  =  -1.593x10”® 

Jg  =  -0.230x10”® 

=  0.502x10”® 

The  values  of  the  geopotential  coefficients  are  those 
derived  by  Kozai  (12) . 

Having  expressed  the  perturbing  acceleration  in  the  STW 
coordinate  frame  the  most  convenient  set  of  perturbation 
equations  is  the  Gaussian  form  of  the  Lagrange  planetary 
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perturbation  equations: 

a  =  (2/n) (l-e^)^/^[esin(f) S+(p/r)T] 
e  =  (1-e^) (na) [sin (f ) S+ (cos (f ) 

+ (e+cos (f ) ) / (1+ecos (f ) ) ) T] 
i  =  rcos(u)/{na^(l-e^)^^^)W 

=  rsin  (u)  /  (na^  (1-e^)  ^'^^sin ( i) )  W  2.16 

oj  =  (1-e^)  V  (nae)  [-cos  (f )  S 
+sin(f) (1+r/p) T] -cos (i)d 
M  =  n+1/ (na^e) [ (pcos (f ) -2er) S- (p+r) sin (f ) T] 


where  S,  T,  and  W  are  given  by  eqs.  2.13-2.15. 

Finally,  eqs.  2.16  are  combined  with  eqs.  2.6  to  form  the 
system  equations  of  motion  (eq.  2.4) . 

In  addition  to  the  zonal  terms  in  the  geopotential 
which  vary  with  latitude,  longitude  dependent  terms, 
called  tesseral  terms,  in  the  geopotential  can  be  included 
in  the  equations  of  motion.  Kaula  (13)  derived  an  expres¬ 
sion  for  the  geopotential,  as  a  function  of  the  classi¬ 

cal  orbital  elements: 


OP 

(i) 


(e)S, 


p^Q^)lmp'*'q=-«''ilpq""“)lmpq 


(oj,M,n  ,9)  , 
2.17 


^here: 


the  F.  (i)  are  the  inclination  functions, 
^mp 

the  Gfl  (e)  are  the  eccentricity  functions, 
ilpq 
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'ilmpq  =  IdT  =os  [  (A-2p)co+ (£-2p+q)  m 

+m(.l-9)]  +  sin([(Z-2p)co 

xm 


+  (Jl-2p+q)M+m(n-e)  ]  ,  and 


the  and  the  are  the  geopotential  coeffi¬ 
cients. 

The  inclination  functions,  eccentricity  functions,  and  the 
geopotential  coefficients  are  given  in  Appendix  A. 

The  variational  equations  of  motion  due  to  the  tesseral 
terms  in  the  geopotential  can  be  calculated  by  applying 
Lagrange's  planetary  equations  due  to  a  disturbing  poten¬ 
tial: 


a  = 


e  = 


1  = 


(2/na)  jjM  , 

(1/na^e)  "T^l  ' 

dn 

3V  3V 

l/(na(l-e^)^/^)  [cot(i)-^  -  csc(i)-^]  , 


csc(i)/{na^(l-e^)  , 

9i 


2.18 


u  =  (l-e^)^/V(na^e) 


3v 


Urn 


3e 


3  V 

-cot(i)/(na^(l-e^)  ,  and 


3i 


3V„ 


M  = 


-(l-e^)/(na^e)  -2/(na)  “ 


3e 


3a 


The  tesseral  harmonics  are  of  small  amplitude  and  to 
first  order  cause  no  secular  variations  in  the  state, 
causing  only  periodic  variations  in  the  classical  orbital 
elements  (10) . 


Two  assumptions  were  made  when  determining  the 
variations  due  to  the  tesseral  harmonics.  First,  only 
the  terms  through  were  included.  The  indices  for  the 
potential,  V,  are  therefore  £=2,3,4  and  m=l,  ...,  £. 

The  terms  in  the  geopotential  for  m=0  are  the  zonal  varia¬ 
tions  and  are  included  in  equation  2.16.  The  second 
assumption  is  that  the  eccentricity  is  small  and  can  be 
neglected  in  the  variational  equations  for  the  tesseral 
harmonic  s . 

Finally,  equations  2.18  are  combined  with  equations 
2.16  and  2.6  to  form  the  system  equations  of  motion 
(eq.  2.4) . 

Other  forces  acting  on  the  satellite,  such  as  solar 
radiation  or  third-body  gravitation  due  to  the  sun  or  the 
moon,  are  not  included  in  the  dynamics  model.  Errors  caused 
by  neglecting  these  forces,  as  well  as  errors  in  the 
dynamics  model  due  to  uncertainties  in  parameters  such 
as  atmospheric  density  and  geopotential  coefficients  will 
be  compensated  for  in  the  process  noise  matrix  described 
in  Chapter  3 . 

2.3.  Observations 

The  observations  used  to  update  the  estimate  of  the 
state  vector  are  slant  range,  azimuth,  and  elevation 
relative  to  a  radar  site.  A  mathematical  model  must  be 
provided  which  relates  the  observations,  z,  to  the  states. 


The  approach  used  in  this  work  is  to  first  express  the 
satellite  position  vector,  r,  in  terms  of  the  state 
vector,  X,  in  the  perifocal  (PQW)  coordin  ::e  system  (14): 


r  =  rcos (f ) P+rsin(f ) Q  2.20 

where  P  and  Q  are  unit  vectors  in  the  orbital  plane;  P  is 
in  the  direction  of  perigee  and  Q  is  perpendicular  to  P  in 
the  direction  of  motion. 

The  next  step  is  to  rotate  the  position  vector  from 
the  perifocal  frame  to  the  geocentric  equatorial  coordinate 
(GEC)  frame  and  then  to  the  topocentric  horizon  coordinate 
(THC)  frame  (14) ; 


^-^THC 


GEC 


“  “2«3 PQW 
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where: 

M^d,!)  =  cos  (fl)  cos  («) -sin  (D)  sin  (cj)  cos  ( i) 
M2(1,2)  =  -cos  (Q)  sin  (oo) -sin  (12)  cos  (oi)  cos  (i) 
M2(1,3)  =  sin(n)sin(i) 

M2(2,1)  =  sin  ((2)  cos  (co) +COS  (12)  sin  (cj)  cos  (i) 
M2(2,2)  =  -sin(n)  sin(co)+cos(C2)cos(oj)cos(i) 


M^O,!)  =  sin(o))  sin(i) 
M2(3,2)  =  cos  (w)  sin  ( i) 
(3,3)  =  cos (i) 


M2  =  -sin(9) 


sin(L)cos(0)  sin(L)sin(9)  -cos (L) 


cos  (9) 


cos(L)cos(9)  cos(L)sin(9)  sin(L) 


L  is  the  geocentric  latitude  of  the  radar  site. 


9  is  the  local  sidereal  time  of  the  radar  site. 


With  the  position  vector,  r,  expressed  in  the  THC 


frame,  the  slant  range  vector  is  now: 

p=r-R=pgS+pgE+p2Z 


2.22 


where  R  is  the  position  vector  of  the  radar  site  on  an 
oblate  Earth  in  the  THC  frame.  Finally,  the  observations, 
z,  can  be  expressed  as  functions  of  the  slant  range  vector. 


t  ^  2  ^  2,1/2 

P  =  (Ps  Pe  Pz^ 

Az  =  sin“^(Pj/(Pg  + 

El  =  sin  ^(pg/p) 


2.23 
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CHAPTER  3 
ESTIMATION 


3.1.  Introduction 

The  primary  purpose  of  estimation  is  to  provide  an 
estimate  of  the  state  vector,  x{t) ,  based  on  a  system 
dynamics  model,  £(x(t) ,t) ,  and  observations  known  to  be 
corrupted  by  noise,  £(t) .  For  the  orbit  determination 
problem  a  continuous,  nonlinear,  systems  model  with  errors 
modeled  by  white-Gaussain  noise,  w(t) ,  is  used  with  dis¬ 
crete  measurements  corrupted  by  white-Gaussian  noise,  v(t) . 
The  system  equations  are; 

x(t)  =  f  (x(t)  , t) +w(t)  ,  w(t)  =  N(0,Q),  3.1 

z(t^)  *  )  ■'■Ik 

v(tj^)  =  V =  N(^R)  .  3.2 

where : 

Q  is  a  diagonal  process  noise  matrix,  and 
R  is  a  diagonal  measurement  noise  matrix. 

Also,  V  and  w  are  assumed  to  be  independent. 

One  method  of  updating  the  state  estimate  with  noisy 
observations  is  a  batch  least  squares  filter.  This  approach 
determines  an  estimate  of  the  state  that  minimizes  the  sum 
of  the  squares  of  the  residuals.  A  residual  is  defined  to 


be  the  value  of  the  actual  observation  minus  the  predicted 
value  based  on  the  state  estimate.  The  batch,  least 
squares  filter  will  provide  the  base-line  filter  for 
evaluating  the  performance  of  other  filters  considered  in 
this  research. 

Two  disadvantages  of  the  least  squares  filter  are 
first  that  it  does  not  use  information  about  the  process 
noise  to  update  the  state.  Second,  it  is  a  batch  processor. 
This  requires  that  sufficient  data  for  an  update  must  be 
stored  prior  to  the  update.  Also,  the  current  estimate 
of  the  state  vector  from  the  batch  processor  does  not  con¬ 
tain  information  based  on  the  most  recent  data,  but  is  based 
on  the  data  prior  to  the  last  update. 

A  second  class  of  filters  considered  in  this  research, 
called  recursive  filters,  include  process  noise  in  their 
state  estimates  and  provide  a  state  estimate  based  on  all  of 
the  data,  up  to  and  including  the  most  recent  observation. 
For  a  linear  problem  with  zero-mean,  white,  Gaussian  noise 
the  optimal  estimator  is  the  Kalman  filter  (15) .  For  the 
nonlinear  problem,  the  Kalman  filter  may  no  longer  provide 
the  optimal  estimate  of  the  state;  however,  the  extended 
Kalman  filter  has  been  successfully  applied  to  several 
nonlinear  problems  described  by  Maybeck  (16) ,  Gelb  (17) , 
and  Jazwinski  (18)  and  will  therefore  be  considered  in  this 
research.  The  extended  Kalman  filter  is  based  on  a  first 


order  Taylor  series  expansion  of  equations  3-1  and  3.2 
about  the  state  estimate.  Retention  of  second  order  terms 
in  the  expansions  results  in  truncated  and  Gaussian  forms 
of  second  order  filters.  Only  the  Gaussian  form  of  the 
second  order  filter  will  be  considered  in  this  research. 

The  following  sections  provide  a  description  of  each 
of  the  filters  investigated  along  with  the  algorithm  used 
for  updating  each  filter. 

3.2.  Least  Squares  Filter 

The  basic  principle  of  the  method  of  least  squares  is 
that  the  best  estimate  of  the  state  vector  is  the  estimate 
which  minimizes  the  svun  of  the  squares  of  the  residuals. 
Define  the  residual  to  be 


3.3 


where : 


2  are  the  actual  observations. 

-o 

Since  the  actual  value  of  the  state,  x,  is  not  known,  h(x) 
must  be  approximated.  This  is  accomplished  using  a  Taylor 
series  expansion  of  h(x)  about  the  nominal  trajectory,  x: 


h(x)  =  b(x)  +  1 (x-&)  +  H.O.T. 


»  z^  +  H6x, 


where; 


3  .4 

3.5 


^  are  the  predicted  values  of  the  observations. 


The  residuals  are  expressed  as: 


r  =  (^-2^)  -H5x  .  3.6 

The  weighted  least  squares  criteria  is  now 

^  (r^Wr)  =  ^  (  [z  -z  -H6x]'^W[z  -z  -H6x]  )=0  3.7 

where 

W  is  the  weighting  matrix. 

The  update  equations  for  the  state  estimate  are 


5x  =  (h'^WH)"^h'^W(z  -z  ) 

—  — o  — c 

3.8 

Vl  " 

3.9 

Because  of  the  nonlinearities  in  the  problem,  this  process 
(eqs  3. 8- 3. 9)  must  be  iterated  until  the  solution  converges, 
i.e.  ^  approaches  zero. 

The  matrix  W  represents  the  confidence,  or  expected 
accuracy,  of  each  observation  relative  to  the  other  obser¬ 
vations.  In  this  study  the  same  weighting  matrix  is  used 
for  all  the  radar  sites  providing  observations  of  slant 
range,  azimuth,  and  elevation.  This  means  that  all  of  the 
slant  range  measurements  are  weighted  equally,  as  well  as 
all  of  the  azimuth  measurements  and  the  elevation  measurements. 
The  weights  used  for  the  slant  range,  azimuth,  and  elevation 
were  based  on  ADCOM  DOS  Technical  Note  79-4  (19)  describing 
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sensor  accuracies.  The  average  standard  deviation  of  the 

error  in  the  range  measurements  was  47  m  or  about  7.4X10~^ 

DU  (1  DU=637a.l45  km).  The  average  standard  deviation  of 

the  error  in  the  azimuth  and  elevation  measurements  was 

-4 

0.014  deg,  or  about  2.4X10  rad.  This  indicates  that  a 
weighting  factor  for  the  slant  range  measurements  should  be 
about  1000  times  the  weighting  factor  of  the  azimuth  and 
elevation  measurements;  however,  faster  convergence  was 
achieved  using  a  ratio  of  100. 

The  H-matrix  is  a  3x7  matrix  relating  changes  in  the 
slant  range,  azimuth,  and  elevation  at  the  time  of  the  ob¬ 
servation  to  changes  in  the  state  at  the  epoch  time,  t^. 
These  partial  derivatives  are  calculated  analytically  rather 
than  numerically  to  reduce  the  computational  load  of  the 
filter.  H  is  evaluated  using  the  chain  rule,  i.e. 


H  =  M^M2M3iyi4M5 


3.10 


where: 


is, .a  3x3  matrix  containing  the  partial  derivatives 
of  slant  range,  azimuth,  and  elevation  with 
respect  to  components  of  a  vector  in  the  THC 


frame  (3) . 

-cos (Az) cos (El)  sin (Az) cos (El) 


sin (El) 


=  sin  (Az)  /  (pcos  (El)  )  cos  (Az)  /  ( pcos  (El)  )  0 

cos (Az) sin (El) /p  -sin (Az) sin (El) /p  cos(El)/p 

3.11 
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M2  is  a  3x3  coordinate  transformation  from  the  THC 
frame  to  the  GEC  frame  (eq.  2.21) 

M2  is  another  coordinate  transformation,  this  time  from 
the  GEC  frame  to  a  coordinate  frame  (UVW)  with 
the  d  and  ^  basis  vectors  in  the  orbital  plane 
and  13  aligned  with  the  position  vector  of  the 
satellite.  M2  is  given  in  equation  2.18  with 
u=f+a)  replacing  o). 

M^  is  a  3x7  matrix  relating  changes  in  the  position 
vector  in  the  UVW  frame  to  changes  in  the  state 
at  some  time,  t.  Since  the  UVW  frame  is  not 
inertial,  changes  in  the  direction  of  the  unit 
vectors  must  be  included.  The  non-zero  elements 
of  M^  are  given  below. 

M^ (1, 1)  =  r/a 
M4(1,2)  =  (r-a)/e 


M4(1,6)  =  aesin (f ) / (1-e^) ^^^-v/n 

M4  (2, 4)  =  rcos  (i) 

M4(2,5)  =  r 

M4(2,6)  =  a^  (l-e^)^^^/r 
M^  (3,3)  =  rsin(u) 


3  .12 


M^(3,4)  =  -rsin ( i) cos (u) 

Mg  is  the  7x7  state  transition  matrix, 
mately  letting 


X  =  +  (t-t  )3t 

—  — o  o  — 


It  is  approxi- 


3.13 


and  evaluating  x  using  values  of  the  state  at  epoch, 
t^ .  Only  secular  and  long-period  variations  in 
the  elements  due  to  drag  and  the  zonal  harmonics 
in  the  geopotential  are  included  in  the  calcula¬ 
tion  of  the  state  transition  matrix. 


Equation  3.13  is  adequate  for  the  variation  in  the 
elements  due  to  the  geopotential,  but  does  not  apply  to  the 
variations  in  the  elements  due  to  atmospheric  drag  due  to 
the  presence  of  position  and  velocity  in  the  variational 
equations.  King-Hele  (8)  shows  that  for  small  (t-t^) , 
the  variations  of  the  elements  due  to  atmospheric  drag  are 
governed  by: 


a  =  a^-QnBa  Ppexp(-z) {I^ (z) +2el^ (z) ) (t-t^) 


e  = 


=  e^-QnBap  exp(-z) (I, (z) -e/2 (I^-l2) ) (t-t^) 


i  =  i^-l/4aa)^BQ^^^p„exp  (-Z)  sin  ( i)  [I  -2el. 
o  e  P  o 


3.14 


0.  = 


0) 


+  (l2-2eI^)cos  (2cu)  ]  (t-t^) 

n^-l/4a(jjgBQ^^^Ppexp  (-Z)  [I^-2eI^]  sin  (  2(o)  (t-t^) 
>1/2 


=  oj^+l/4a«gBQ  ppexp(-z)cos(i)  [I^-2eI^] 


sin ( 2«) ( t-t  ) , 
o 


where: 


p  is  the  density  at  perigee, 
z  =  ae/H,  and 

the  are  the  modified  Bessel  Functions  of  order 
The  state  transition  matrix  is  given  in  Appendix  B. 

The  algorithm  used  to  estimate  the  state  using  a  least 
squares  filter  is  presented  below. 

1.  Select  a  batch  size,  m,  for  processing  the  data. 
Batch  sizes  of  15  -  25  sets  of  observations  (slant 
range,  azimuth,  and  elevation)  are  used  in  this 
research. 

2.  Read  in  the  next  m  sets  of  observations  and  save 

as  z^ . 

-o 

3.  Propagate  the  state,  x,  to  each  observation  time, 
t . ,  using  eq  2.4  and  calculate  z  (t.)  and  H{t.). 

•  “O  X  X 

Combine  the  ^(t.)  's  and  H(t.)  's  to  form  z  (t) 
and  H  ( t)  . 

4.  Calculate  ^  from  eq  3.8. 

5.  Update  x,  eq.  3.9. 

6.  Repeat  3-5  until  solution  converges. 

7.  Propagate  x,  eq  2.4,  to  new  epoch  time  (time  of 
last  processed  observation)  and  read  in  next  m 
data  sets. 

8.  Repeat  until  all  of  the  data  are  processed. 


The  Kalman  filter  is  a  sequential  estimator  providing 
an  estimate  of  the  state  vector  based  on  the  system  dynamics 
model  and  information  from  all  of  the  previous  observations. 
For  a  linear  problem  with  Gaussian  additive  noise  the 
Kalman  filter  provides  an  optimal  estimate  of  the  state 
vector.  However,  for  the  nonlinear  problem,  the  Kalman 
Filter  no  longer  guarantees  an  optimal  solution.  This  is 
due  to  the  fact  that  the  optimal  estimate  requires  knowledge 
of  the  probability  density  function  of  the  states,  p(x,t). 
For  a  nonlinear  problem,  p(x,t)  will  not  remain  a  Gaussian 
density  function  even  though  the  initial  state  may  be 
Gaussian  and  the  process  noise  and  observation  noise  are 
Gaussian.  If  £(x,t)  is  not  Gaussian,  the  estimation  of  the 
optimal  state  vector  requires  the  knowledge  of  the  entire 
density  function,  which  is  not  feasible.  Therefore,  some 
form  of  a  suboptimal  estimator  is  required.  The  simplest 
form  is  the  linearized  Kalman  filter.  The  linearized  Kalman 
Filter  (16)  provides  an  optimal  estimate  of  the  error  in 
the  state,  ^(t).  This  error  is  then  added  to  the  nominal 
state  x(t)  to  determine  the  total  state: 

x(t)=x  (t)+  ^(t)  3.15 

This  form  of  the  Kalman  filter  is  adequate  if  the  nominal 
state  is  nearly  the  same  as  the  "true"  state;  however,  if 
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the  nominal  state  and  the  "true"  state  differ  by  a  large 
amount  the  linearization  is  no  longer  appropriate. 

The  extended  Kalman  filter  is  an  attempt  to  reduce 
the  effect  of  this  problem.  The  EKF  relinearizes  about 
the  "new"  state  after  each  measurement  update.  This  reduces 
the  difference  between  the  optimal  state  and  the  "true"  state 
thus  enhancing  the  assumption  that  a  linearization  of  the 
dynamics  and  observations  is  valid.  The  linearized  system 
dynamics  and  observation  equations  are  obtained  by  using  a 
Taylor  series  expansion  of  the  equation  about  the  nominal 
trajectory.  A  derivation  of  the  EKF  equations  is  presented 
by  Maybeck  (16)  and  will  not  be  presented  here.  The  result¬ 
ing  EKF  equations  are  (16) : 

Propagation: 

k(t)  =  f  (&(t)  ,t)  , 

P(t)  =  F(x,t)P(t)  +  P(t)F'^(&,t)  +  Q,  3.16 

Update : 

K  =  P(-)H'^(k(-)  )  (H(k(-)P(-)H'^(k(-)  )  +  R]~^, 

ig(  +  )  =  k(-)  +  K[z  -  h(k(-))],  and 
—  —  — o  —  — 

P(  +  )  =  P(-)  -  KH(k(-)  )  P(-)  ,  3.17 

where: 


P  is  the  state  covariance  matrix, 

F  is  the  linearized  dynamics  matrix, 
Q  is  the  process  noise  matrix, 

K  is  the  filter  gain  matrix. 


: 


r 

‘5 
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a 
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H  is  the  linearized  observation  matrix,  and 


R  is  the  measurement  noise  matrix. 


The  initial  values  of  the  state  and  covariance  matrix 


must  be  input  to  the  filter.  The  initial  values  of  the 


state,  called  the  epoch  values,  are  described  in  Chapter  4 


for  each  case.  The  initial  value  of  the  covariance  is  a 


measure  of  the  confidence  in  the  epoch  values  of  the  state. 


The  initial  covariance  matrix  is  assumed  to  be  a  diagonal 


matrix  with  each  of  the  diagonal  elements  equal  to  lE-06. 


The  process  noise  matrix,  Q,  is  a  measure  of  the  uncertainty 


in  the  dynamics  model.  Sources  of  error  in  the  dynamics 


include  unmodeled  forces  such  as  third  body  perturbations 


or  higher  order  terms  in  the  geopotential,  uncertainties  in 


the  modeled  geopotential  coefficients,  and  errors  in  the 


atmospheric  density  model.  The  process  noise  matrix  used 


in  this  research  is  described  in  Chapter  5.  The  H-matrix 


is  the  same  as  the  H-matrix  used  in  the  least  squares 


filter  and  is  described  in' the  previous  section.  Finally, 


the  measurement  noise  matrix  is  a  measure  of  the  accuracy  of 


the  observations:  slant  range,  azimuth,  and  elevation. 


The  EKF  equations  allow  the  flexibility  to  use  different 


weights  for  each  type  of  observation;  slant  range,  azimuth. 


or  elevation,  and  for  each  radar  site.  For  the  radar  sites 


used  in  this  work,  the  standard  deviations  of  the  measure¬ 


ments  did  not  vary  significantly  from  site  to  site  (19)  and 


therefore  only  one  measurement  noise  matrix  is  used.  R  is 

assumed  to  be  a  diagonal  matrix.  The  diagonal  elements 

are  equal  to  the  variances  in  the  slant  range,  azimuth, 

and  elevation  measurements  respectively.  Initial  values 

for  the  variances  were  calculated  by  squaring  the  average 

of  the  standard  deviations  for  each  of  the  radar  sites. 

The  average  standard  deviation  of  the  slant  range  measure- 

ments  is  47  m  which  results  in  a  variance  of  about  lE-11  DU  . 

Similarly,  the  average  standard  deviations  of  the  azimuth 

and  elevation  are  0.14  deg  which  result  in  variances  of 
2 

5E-08  rad  .  The  values  of  the  elements  of  the  R-matrix 
were  later  altered  to  improve  the  filter  performance.  The 
values  used  are  reported  in  the  results  section.  Chapter  5. 

Several  modifications  to  the  EKF  can  be  incorporated 
to  improve  the  filter  performance.  One  change  is  the 
inclusion  of  "pseudonoise”  (16)  in  the  Q-matrix  and  the 
R-matrix  tO'  compensate  for  unmodeled  errors  and  errors  due 
to  nonlinearities  in  the  problem.  This  method  is  used  to 
alter  the  filter  gain  and  achieve  a  balance  between  ignoring 
the  new  data  (low  filter  gain)  and  tracking  the  noise  in 
the  data  (high  filter  gain).  Another  modification  to  com¬ 
pensate  for  nonlinearities  in  the  problem  is  to  add 
iterations  in  the  propagation  and  update  equations.  An 
iterated,  extended  Kalman  filter  which  incorporates 
iterations  in  the  update  equations  is  investigated  in 
this  research.  The  modified  update  equations  are  (16) : 


where: 


P(-)  is  the  covariance  just  before  the  update, 

and 

x(l)  is  the  estimate  from  eq.  3.17, 

=  x(-)  +  K[z  -  h(x(-)  )]  . 

—  -o  —  — 

Ec[uations  3.18  are  iterated  until  the  solution  converges  or 
a  maximum  number  of  iterations  is  reached.  After  the  itera¬ 
tion  of  equations  3.18,  the  new  P-matrix  is  calculated  from 
equation  3.17. 

The  algorithm  used  for  the  extended  Kalman  filter  is 
listed  below. 

1.  Read  in  the  initial  state  vector,  covariance  matrix 
and  epoch  time. 

2.  Read  in  the  first  data  set. 

3.  Propagate  the  state  vector  and  covariance  matrix  to 
the  time  of  the  first  observation  (eq.  3.16). 

4.  Calculate  the  H-matrix  and  F-matrix. 

5.  Calculate  the  gain  matrix  and  update  the  state 
vector  (eq.  3.17). 

6.  Repeat  4  and  5  until  the  solution  converges  using 
equation  3.18  to  calculate  the  new  gain  and  the 


new  state  vector. 

7.  Calculate  the  updated  covariance  matrix  using 
equation  3.17. 

8.  Update  the  epoch  time  to  the  time  of  the  last 
observation. 

9.  Read  in  next  data  set  and  repeat  steps  3-9  until 
all  of  the  data  are  processed. 

3.4.  Second  Order  Filters 

If  the  nonlinearities  in  the  problem  are  significant, 
inclusion  of  second  order  terms  in  the  Taylor  series 
expansions  of  the  dynamics  and  observation  equations  along 
with  assumptions  about  the  conditional  density  function  of 
the  state  based  on  the  observations  should  improve  the  filter 
performance  (16) .  Two  filters  of  this  type  are  developed  by 
Maybeck  (16) .  Only  the  Gaussian  form  of  the  second  order 
filter  will  be  investigated. 

As  stated  before,  the  optimal  estimate  of  the  state, 

/N 

x(t),  based  on  all  of  the  previous  measurements  requires 
the  knowledge  of  the  entire  conditional  probability  density 
function  of  the  state  based  on  the  measurements.  Since  this 
is  not  feasible,  certain  assumptions  about  the  density 
function  must  be  made.  One  set  of  assumptions  is  to  first 
assume  that  the  density  function  is  symmetric  about  its 
mean,  thus  all  of  the  odd  moments  of  order  three  and  higher 
are  zero.  Also,  assume  that  the  density  function  is 
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concentrated  about  the  mean  so  that  the  fourth  order  and 
higher  order  even  moments  may  be  neglected.  Combining 
these  assumptions  with  the  inclusion  of  the  second  order 
terms  in  the  Taylor  series  expansions  of  the  dynamics  and 
observation  equations  results  in  the  truncated  second 
order  filter.  If  the  additional  assumption  that  the  pro¬ 
cess  noise  is  time  invariant  is  made,  the  resulting 
modified  truncated  second  order  filter  equations  are 
(16)  : 


Propagation: 

x(t)  =  f[x(t),t]  + 


Tr- 


P(t)  =  FCx(t)  ,t]P(t)  +  P(t)F^[x(t)  ,t]  +  Q,  3.19 


where: 


F  and  Q  are  the  same  as  in  the  EKF,  and 


^(t)  is  an  n-vector  whose  kth  element  is: 


d^f, 


bdk(t)  =  (1/2)  tr{ 


— ^  I  P(t)}; 

dx^  x=k 


Update : 


A(t)  =  H[x(t-)  ,tlP(t-)H  [x(t-)  ,t] 


-b  (t-)^(t-)  +  R(t)  , 
— m  — m 


K(t)  =  P(t-)H’^(x(t-)  ,t]A"^(t)  , 
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where: 


x(t+)  =  x(t-)  +  K(t)  [z  -z  -b  (t-)  ]  ,  and 
—  —  -o  — c  — m 

P(t+)  =  P(t-)  -  K(t)H[x(t-)  ,  t]P{t-)  , 


H  and  R  are  the  same  as  in  the  EKF,  and 

b  (t)  is  an  m-vector  whose  kth  element  is; 
— m 
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3x  x=x 


P(t-)  }. 


Similarly,  the  Gaussian  second  order  filter  is  derived 
by  assuming  that  the  conditional  probability  density  function 
remains  Gaussian.  The  third  order  and  higher  order  odd 
moments  of  a  Gaussian  density  function  are  zero  and  the 
fourth  order  and  higher  order  even  moments  can  be  expressed 
in  terms  of  the  covariance.  If  only  the  fourth  order  moment 
is  considered  and  the  process  noise  matrix  is  assumed  to 
be  time  invariant,  the  resulting  modified  Gaussian  second 
order  filter  equations  are  the  same  as  the  truncated  second 
order  filter  with  the  exception  of  the  filter  gain,  A(t) , 
which  is  now  (16)  : 


where: 


A(t)  =  H[x(t-)  ,tlP(t-)  [x(t-)  ,  t] 
+  B(t-)  +  R(t)  , 


B(t-)  is  an  m-by-m  matrix  with  kith  element: 


3.21 


Bki(t-)  =  (1/2) tr{ 


P(t-) 


P(t-)} . 


3x  x=sx 


3x  x=x 


,  t.  > 


The  differences  between  the  EKF  and  the  second  order 
filters  are  the  additional  terms  in  the  dynamics  equations, 
state  update  equations,  and  the  filter  gain.  These  addi¬ 
tional  terms  in  the  filter  equations  can  be  considered  bias 
terms.  These  bias  terms  can  be  included  individually  into 
the  filters,  making  it  possible  to  investigate  the  effect 
of  each  term  separately.  The  dynamics  and  update  bias 
terms  are  the  same  for  the  two  second  order  filters,  so  the 
only  difference  between  the  two  second  order  filters  is  the 
bias  term  in  the  filter  gains. 

The  calculation  of  the  bias  terms  in  the  second  order 

filters  requires  the  evaluation  of  the  second  partial 

derivatives  of  the  states  at  time  t,  x(t) ,  with  respect 

to  the  states  at  time  t  ,  x(t  ) .  Since  these  terms  are 

0  —  0 

already  second  order,  only  the  two-body  dynamics  plus 
secular  variations  due  to  the  Earth's  oblateness,  J2»  and 
atmospheric  drag  are  considered  in  the  bias  terms.  The 
derivations  of  the  second  order  bias  terms  used  in  this 
research  are  presented  in  Appendix  C. 

The  algorithm  used  for  the  second  order  filters  is 
the  same  as  the  algorithm  used  for  the  EKF  with  the  inclu¬ 
sion  of  the  bias  terms  where  appropriate. 

The  system  dynamics  and  observation  equations  described 
in  Chapter  2  can  now  be  combined  with  the  filter  equations 
described  in  this  chapter  to  solve  the  orbit  estimation 
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problem.  The  four  cases  used  to  compare  the  performance 
of  the  various  filters  are  presented  in  the  next  chapter, 


CHAPTER  4 


OBSERVATIONAL  DATA 

4.1.  Introduction 

The  objective  of  this  research  is  to  compare  the  per¬ 
formance  of  several  filters  using  actual  satellite  tracking 
data.  The  first  step  in  the  process  is  to  verify  the  filter 
algorithms.  This  is  accomplished  using  synthetic  data 
which  is  free  of  process  noise  and  measurement  noise.  In 
addition  to  this  test  case,  observational  data  obtained 
from  tracking  three  different  satellites  are  processed  in 
order  to  compare  the  performance  of  the  filters.  The  first 
satellite  orbit  is  a  high  altitude  (1000  km  at  perigee) , 
non-circular  (e=0.015),  polar  (i=90  deg)  orbit.  The 
second  orbit  is  a  low  altitude  (2  50  Ion  at  perigee)  ,  non¬ 
circular  (e=0.01),  inclined  (i=72  deg)  orbit.  The  third 
satellite  orbit  is  a  low  altitude  (300  km) ,  nearly  circular 
(e=0.0003)  orbit  with  an  inclination  of  50  degrees. 

The  epoch  values  for  the  test  case  and  the  three  satel¬ 
lite  orbits  are  presented  in  the  next  section  and  a  discus¬ 
sion  of  the  data  is  presented  in  Section  4.3. 

4.2.  Epoch  Values  of  the  Satellite  Orbits 

The  epoch  values  for  the  test  case  and  the  three 
satellite  orbits  considered  in  this  research  are  presented 


38 


in  Tables  4-1  through  4-4.  The  epoch  values  for  satellites 
4507,  6633,  and  1  299  are  the  singularly  averaged  classical 
eleinents  provided  by  the  Air  Force  SPACECOM/DOA .  The  epoch 
values  for  the  test  case  are  based  on  satellite  4507  with 
small  variations  in  the  elements  representing  initial 
errors  in  the  state  estimate. 

Table  4-1 

Epoch  Values  for  the  Test  Case 


Time: 


4:24:28.999  on  day  19,  1984. 


Semi-major  axis: 

Eccentricity: 

Inclination: 

Longitude  of  Ascending 
Node: 

Argument  of  Perigee: 
Mean  Anomaly: 

Ballistic  Coefficient: 


1.1703  DU 
0.0147 
90.076  deg 

293.363  deg 
214.900  deg 
316.080  deg 
0.00698  m^/kg 


Table  4-2 


Epoch  Values  for  Satellite  4507 


Time; 

Semi-major  Axis: 

Eccentricity: 

Inclination: 

Longitude  of  Ascending 
Node: 

Argument  of  Perigee; 
Mean  Anomly: 

Ballistic  Coefficient: 


4:24:28.999  on  day  19,  1984 
1.1703  DU 
0.01475 
90. 077  deg 

293.363  deg 
214.908  deg 
316.078  deg 
0.00698  m^/kg 


Table  4-3 

Epoch  Values  for  Satellite  10299 


Time: 

Semi-major  Axis; 

Eccentric ity : 

Inclination ; 

Longitude  of  Ascending 
Node : 

Argument  of  Perigee 
Mean  Anomaly: 

Ballistic  Coefficient; 


23:41:24.000  day  245,  1977 
1.0401  DU 
0.00955 
72.844  deg 

115.989  deg 
56.617  deg 
106.653  deg 
0.0010  m^/kg 


Table  4-4 


Epoch  Values  for  Satellite  6633 


Time: 

Semi -major  Axis: 

Eccentric ity : 

Inclination: 

Longitude  of  Ascending 
Node : 

Argument  of  Perigee: 
Mean  Anomaly: 

Ballistic  Coefficient: 


22:44:40.000  day  69,  1979 
1.0545  DU 
0.0003 
50.500  deg 

283.827  deg 
169.807  deg 
348.612  deg 
0.0120  mVkg 


4.3.  Data  Summary 

The  observational  data  used  in  this  research  consists 
of  radar  passes  and  observations.  A  radar  pass  is  the  hori¬ 
zon  to  horizon  coverage  of  a  satellite  orbit  by  a  particular 
radar  site.  An  observation  is  the  measured  slant  range, 
azimuth,  and  elevation  of  a  satellite  with  respect  to  a 
particular  radar  site  at  a  given  time.  The  radar  site 
locations,  provided  by  Air  Force  SPACECOM/DOA,  are  given 
in  Appendix  D. 

The  observational  data  for  the  test  case  are  synthetic 
data  based  on  the  epoch  values  of  satellite  4507.  The  data 
are  generated  by  propagating  the  initial  state  forward  using 
the  dynamics  described  in  Section  2.2.  When  the  satellite 
is  visible  to  a  radar  site  a  set  of  look  angles;  predicted 


slant  range,  azimuth,  and  elevation,  are  calculated.  For 
this  synthetic  data,  a  radar  pass  consists  of  ten  sets  of 
look  angles,  each  separated  by  60  seconds.  These  look 
angles  are  saved  and  are  used  as  noise  free  observational 
data.  There  are  a  total  of  100  observations;  ten  radar 
passes  with  ten  observations  per  pass,  over  a  24  hour 
period  in  the  test  data.  This  gives  one  radar  pass  every 
1.35  satellite  revolutions.  The  maximum  time  between 
observations  is  5.18  hours. 

The  data  for  the  other  three  satellites  are  actual 
observational  data  provided  by  Air  Force  SPACECOM/DOA.  The 
data  include  time  of  observation,  satellite  niamber,  radar 
site  number,  and  slant  range,  azimuth,  and  elevation  of 
the  satellite  relative  to  the  site.  The  radar  passes 
vary  from  two  to  15  observations  per  pass  and  the  obser¬ 
vations  in  a  particular  pass  are  separated  by  six  seconds 
up  to  one  minute.  Also,  due  to  the  geometry  of  the  orbits 
and  tasking  requirements,  radar  passes  are  separated  by 
from  a  few  minutes  up  to  seven  hours. 

Three  days  of  data,  days  19  -  21  in  1984,  containing 
35  radar  passes,  are  processed  for  satellite  4507.  There 
are  a  total  of  215  sets  of  observations,  resulting  in  an 
average  of  one  radar  pass  every  1.18  satellite  revolutions 
and  an  average  of  6.1  sets  of  observations  per  radar  pass. 
The  maximum  time  between  observations  is  5.33  hours. 


One  day  of  observational  data,  day  246,  1984,  contain¬ 
ing  15  radar  passes  and  72  sets  of  observations  are  pro¬ 
cessed  for  satellite  10299.  This  results  in  an  average  of 
one  radar  pass  every  1.07  satellite  revolutions  and  an 
average  of  4.8  sets  of  observations  per  radar  pass.  The 
maximum  time  between  observations  is  4.63  hours. 

Finally,  the  observational  data  for  satellite  6633 
consists  of  12  radar  passes  with  68  sets  of  observations 
from  22:44  on  day  69  through  day  70,  1979.  This  results  in 
an  average  of  one  radar  pass  every  1.37  satellite  revolutions 
with  an  average  of  5.7  observations  per  radar  pass.  The 
maximum  time  between  obseirvations  for  satellite  6633  is 
7 . 2  hour  s . 

The  data  are  summarized  in  Table  4-5. 

Table  4-5 

Observational  Data 


Satellite 

Hours 

Radar 

Passes 

Number 
of  Obs 

Revs/ 

Pass 

Obs/ 

Pass 

Maximum 

Outage 

Test 

24 

10 

100 

1.35 

10 

5.18 

4507 

72 

35 

215 

1.18 

6.1 

5.33 

10299 

24 

15 

72 

1.07 

4.8 

4.63 

CHAPTER  5 


RESULTS 


5.1.  Introduction 

The  performance  of  the  filters  investigated  in  this 
research  are  compared  using  four  criteria.  These  criteria 
are  estimation  errors,  prediction  errors,  computer  time 
of  operation,  and  computer  storage  requirements. 

Estimation  errors  are  based  on  the  state  estimate  at 
time  t(i)  using  all  of  the  observations  up  to  and  including 
time  t(i) .  A  measure  of  the  filter  estimation  errors  is 
obtained  by  comparing  the  residuals  (observed  value  minus 
predicted  value)  in  the  observations;  slant  range,  azimuth, 
and  elevation,  for  each  filter.  Estimation  errors  are 
compared  in  Section  5.2. 

Similarly,  prediction  errors  are  based  on  the  state 
estimate  at  time  t(i)  based  on  observations  up  to  the  time 
t(j)  where  t(i)  is  greater  than  t(j) .  Again,  a  measure  of 
the  predction  errors  is  obtained  by  comparing  the  observation 
residuals  of  the  filters.  Prediction  errors  are  compared 
in  Section  5.3. 

Finally,  computer  time  of  operation  and  computer  storage 
requirements  are  discussed  in  Sections  5.4  and  5.5  respec¬ 
tively. 


5.2.  Estimation  Errors 

The  estimation  errors  of  the  filters  are  compared  by 
comparing  the  residuals  in  slant  range,  azimuth,  and  eleva¬ 
tion  for  each  filter.  The  filters  all  start  with  the  same 
epoch  values  for  each  satellite  and  process  three  days  of 
data  for  satellite  4507  and  one  day  of  data  for  the  test  case 
(synthetic  data  based  on  orbit  of  satellite  4507),  satellite 
10299,  and  satellite  6633.  As  the  data  are  processed  the 
value  of  the  magnitude  of  the  residuals  based  on  the  state 
estimate  after  the  final  iteration  are  saved  and  used  to 
calculate  the  average  residuals  for  each  filter.  Only  the 
last  observation  residual  from  each  radar  pass  is  used  to 
calculate  the  average  residual. 

The  accuracy  and  stability  of  the  Kalman  filter  and 
Gauss  filter  depend  upon  the  values  of  the  initial  covari¬ 
ance  matrix,  process  noise  matrix,  and  measurement  noise 
matrix.  If  these  matrices  are  not  modeled  correctly  two 
problems  can  occur.  One  possibility  is  that  the  estimates 
of  the  standard  deviations  of  the  state  vector,  the  diagonal 
elements  of  the  covariance  matrix,  get  very  small.  This 
implies  that  the  filter  "knows”  the  state  very  well  and 
ignores  the  new  data.  The  true  standard  deviations  of  the 
state  vector  become  much  larger  than  the  estimated  standard 
deviations  and  the  filter  subsequently  diverges.  The 
opposite  situation  occurs  when  the  estimated  standard 


deviations  of  the  state  are  much  larger  than  the  true 
standard  deviations  of  the  state.  In  this  case  the  filter 
estimate  begins  tracking  the  noise  in  the  data  and  the  filter 
overcompensates  for  the  errors  in  the  state  estimate.  Values 
of  the  initial  covariance  matrix,  process  noise  matrix,  and 
measurement  noise  matrix  are  adjusted  by  trial  and  error 
until  satisfactory  filter  performance  is  achieved. 

The  initial  value  of  the  covariance  matrix  is  based  on 
the  confidence  of  the  initial  state  vector.  The  initial 
state  vector  is  the  singularly  averaged  classical  orbital 
element  set  provided  by  SPACECOM/IDOA.  The  state  vector  in 
this  research  consists  of  the  osculating  classical  orbital 
elements  which  differ  from  the  singularly  averaged  elements 
by  the  inclusion  of  short  period  variations  in  the  osculating 
elements.  Therefore,  the  errors  in  the  initial  elements  are 
due  to  differences  between  these  two  sets  of  elements.  For 
example,  the  magnitude  of  the  short  period  variation  in 
semi-major  axis  due  to  J2  is  about  8  km,  or  l.OE-03  DU. 

This  results  in  a  standard  deviation  of  l.OE-06  for  the 
semi-major  axis.  Similar  calculations  for  the  other  ele¬ 
ments  provide  the  remaining  diagonal  elements  for  the  initial 
covariance  matrix.  Starting  with  these  values  the  diagonal 
elements  of  the  initial  covariance  matrix  are  then  varied 
by  trial  and  error  until  satisfactory  filter  performance 
is  achieved.  The  values  for  the  initial  covariance  matrix 


for  the  Kalman  filter  and  Gauss  filter  used  for  the  follow¬ 
ing  comparisons  are: 

P  =  diag(1.0E-06,. . .,1.0E-06)  5.1 

2  2 
where  the  units  are  DU  for  the  semi-major  axis,  rad  for 

2  2 

the  angles,  and  (m  /kg)  for  ballistic  coefficient. 

The  process  noise  matrix,  Q,  is  a  function  of  the 
unmodeled  dynamics  and  uncertainties  in  the  parameters  used 
in  the  dynamics  model.  The  initial  values  of  the  process 
noise  matrix  are  based  on  the  values  used  by  Taylor  (6) . 

The  values  are  again  adjusted  by  trial  and  error  until 
satisfactory  performance  is  achieved.  The  values  used  in 
the  filters  are; 

Q  =  diag(1.0E-10,lE-16, . . . ,lE-06) .  5.2 

2 

The  units  for  the  elements  of  Q  are  (DU)  /TU  for  the  semi- 

2 

major  axis,  1/TU  for  eccentricity,  (rad)  /TU,  for  the  angles, 
2  2 

and  (m  /kg)  /TU  for  the  ballistic  coefficient. 

The  measurement  noise  matrix,  R,  represents  the  uncer¬ 
tainty  in  the  observations.  Starting  values  for  R  are  based 
on  sensor  statistics  described  in  Chapter  2.  The  same  matrix 
is  used  for  all  of  the  radar  sites.  The  measurement  noise 
matrix  used  here  is; 

R  =  diag(1.0E-10,1.0E-07,1.0E-07) .  5.3 

2  2  2 
The  units  are  (DU)  ,  (rad)  ,  and  (rad)  . 


The  first  satellite  data  investigated  in  this  research 
is  the  test  case  used  to  verify  the  filter  algorithms  which 
contains  no  process  noise  and  no  measurement  noise.  The 
filter  estimation  errors  are  presented  in  Table  5-1. 


Table  5-1 

Estimation  Residuals  of  the  Test  Case 


Filter 

SR  Res  (m) 

Az  Res (deg) 

El  Res (deg) 

LSF 

129 

0.012 

0.0047 

EKF 

96 

0.0028 

0.0028 

GSF 

96 

0.0027 

0.0029 

The  errors  in  Table  5-1  are  due  in  part  to  the  evalua¬ 
tion  of  the  H-matrix  and  the  time  between  observations.  In 
the  formulation  of  the  H-matrix,  it  was  assumed  that  the 
orbital  elements  were  nearly  constant  over  the  period 
between  the  observations  so  that 

x(t)  =  x(t  )  +  x(t  ) (t-t  ) .  5.4 

—  —  o  —  o  o 

This  approximation  is  valid  for  small  (t-t^)  or  if  the  mag¬ 
nitudes  of  the  variation  of  the  elements,  x(t) ,  are  small. 
However,  the  average  time  between  observations  is  a  little 
more  than  two  hours,  or  1.35  satellite  revolutions.  Over 
this  time  period,  the  elements  propagate  through  more  than 
one  complete  revolution  of  the  short  period  variations. 

The  magnitude  of  the  short  period  variations  are  large 
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enough  to  reduce  the  accuracy  of  the  filters.  Since  there 
is  no  process  noise  and  no  measurement  noise  in  the  data  for 
the  test  case,  the  errors  in  Table  5-1  are  representative 
of  the  errors  caused  by  the  assumptions  made  in  the  use 
of  equation  5-4. 

The  next  satellite  investigated  in  this  research  is 
satellite  4507.  The  orbit  of  satellite  4507  is  nearly  polar, 
has  an  eccentricity  of  0.015,  and  an  altitude  of  1000  km. 

At  this  altitude  atmospheric  drag  has  very  little  effect  on 
the  orbit.  The  residuals  for  satellite  4507  are  presented 
in  Table  5-2. 


Table  5-2 

Estimation  Residuals  for  Satellite  4507 


Filter 

SR  Res  (m) 

Az  Res  (deg) 

El  Res 

LSF 

140 

0.012 

0.014 

EKF 

425 

0.070 

0.029 

GSF 

426 

0.07  0 

0.028 

The  residuals  in  Table  5-2  are  based  on  three  days  of 
observations  and  are  the  mean  of  the  magnitudes  of  the  last 
observation  residual  for  each  radar  pass. 

The  difference  in  the  magnitude  of  the  residuals 
between  the  least  squares  filter  and  the  recursive  filters 
is  due  in  part  to  the  problem  of  trying  to  estimate  an 
unobservable  state.  At  an  altitude  of  1000  km  the  atmospheric 
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density  is  nearly  zero  so  that  the  effect  of  atmospheric 
drag  is  insignificant.  The  orbit  of  the  satellite  is 
independent  of  the  ballistic  coefficient  so  the  observa¬ 
tions  of  the  orbit  provide  no  information  about  the  ballis¬ 
tic  coefficient.  This  can  be  checked  numerically  by  check¬ 
ing  the  condition  number  of  the  covariance  matrix  (20) . 

The  condition  number  of  a  square,  symmetric  matrix  is  a 
scalar  which  is  the  ratio  of  the  maximum  eigenvalue  to  the 
minimum  eigenvalue  of  the  matrix.  The  magnitude  of  the 
condition  number  is  a  measure  of  the  numerical  stability  of 
the  algorithm;  the  larger  the  condition  number  the  greater 
the  loss  of  accuracy  in  the  estimation  problem.  The  condi¬ 
tion  number  of  the  covariance  matrix  for  the  Extended 
Kalman  filter  for  satellite  4507  is  approximately  lE+10. 

It  is  difficult  to  say  much  about  this  case  alone,  but  the 
condition  number  for  this  case  can  be  compared  to  the  condi¬ 
tion  number  for  satellite  10299,  where  the  ballistic  coeffi¬ 
cient  should  be  more  observable.  The  condition  number  of 
the  covariance  matrix  for  satellite  10299  is  lE+06. 

The  condition  number  itself  does  not  provide  information 
about  which  state  is  unobservable;  however,  for  satellite 
4507  it  is  reasonable  to  assume  that  the  unobservable  state 
is  the  ballistic  coefficient.  It  is  also  reasonable  to 
reduce  the  order  of  the  filter  in  this  case  since  assuming 
that  the  ballistic  coefficient  is  a  constant  and  not 
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estimating  it  will  not  significantly  alter  the  system 
dynamics;  however,  it  should  remove  the  unobservable  state 
from  the  problem.  The  orders  of  both  the  Kalman  filter 
and  the  Gauss  filter  were  reduced  by  not  estimating  the 
ballistic  coefficient  for  satellite  4507.  The  results  of 
the  estimation  problem  are  presented  in  Table  5-3. 


Table  5-3 


Estimation  Residuals  for  Satellite  4507 
Using  Reduced  Order  Filters 


V 

X 

Filter 

SR  Res 

(km) 

Az  Res  (deg) 

El  Res 

EKF 

341 

0.069 

0.03  0 

1 

GSF 

341 

0.069 

0.030 

I 

The  condition 

number 

for 

the  reduced  order 

Kalman 

covariance  matrix  is  lE+07. 

The  next  satellite  orbit  investigated  is  that  of 


satellite  10299,  in  a  low-altitude  (h^  =  250  km),  non¬ 


circular  (e  =  0.01)  orbit,  with  an  inclination  of  72 
degrees.  The  estimation  residuals  for  satellite  10299  are 
shown  in  Table  5-4. 

Table  5-4 

Estimation  Residuals  for  Satellite  10299 


Filter 

SR  Res  (m) 

Az  Res  (deg) 

El  Res  (deg) 

LSF 

185 

0.048 

0.035 

EKF 

221 

0.059 

0.057 

GSF 

189 

0.052 

0.041 
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The  condition  number  for  satellite  10299  is  lE+06. 

The  improvement  of  the  Gauss  filter  over  the  Kalman 
filter  can  be  further  investigated  by  considering  each  bias 
term  individually.  The  bias  terms  in  the  Gauss  filter 
have  been  described  in  Chapter  3.  There  are  bias  terms  in 
the  dynamics,  observations,  and  filter  gain.  The  data  for 
satellite  10299  is  processed  using  the  Gauss  filter,  first 
with  no  bias  terms,  then  with  only  the  dynamic  bias  term, 
the  observation  bias  term,  and  the  gain  bias  term  individ¬ 
ually.  The  results  are  presented  in  Table  5-5. 

Table  5-5 


Effect  of  Bias  Terms  in  the  Gauss  Filter 


for 

Satellite 

10299 

Bias  Term  SR 

Res  (m) 

Az  Res  (deg) 

El  Res  (deg) 

None 

221 

0.059 

0.057 

Dynamics 

221 

0.059 

0.057 

Observations 

189 

0.052 

0.041 

Gain 

220 

0.059 

0.057 

All 

189 

0.052 

0.041 

The  results  indicate  that  the  bias  term  in  the  observa¬ 
tion  equations  provides  the  most  significant  improvement  for 
the  Gauss  filter.  It  is  also  noted  that  the  Causs  filter 
with  no  bias  terms  reduces  to  the  Kalman  filter. 

The  final  satellite  case  considered  is  satellite  6633, 
which  is  in  a  nearly  circular  {e=0.0003),  low  altitude 
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(h=300  km)  orbit  with  an  inclination  of  50  degrees.  Satel¬ 
lite  6633  is  known  to  be  Skylab  which  is  non-symmetric 
and  tumbling;  therefore,  the  ballistic  coefficient  is  not 
constant.  The  estimation  residuals  for  satellite  6633  are 
presented  in  Table  5-6. 


Table  5-6 

Estimation  Residuals  for  Satellite  6633 


Filter 

SR  Res  (m) 

Az  Res  (deg) 

El  Res  (deg) 

LSF 

106 

0.028 

0.013 

EKF 

512 

0.098 

0.151 

GSF 

474 

0.090 

0.150 

As  was  done  for  satellite  10299,  the  bias  terms  for  the 
Gauss  filter  can  be  investigated  individually  for  satellite 
6633.  The  results  are  presented  in  Table  5-7. 

Table  5-7 

Effect  of  Bias  Terms  in  the  Gauss  Filter 
for  Satellite  6633 


;ias  Term 

SR  Res  (m) 

Az  Res  (deg) 

El  Res  (deg) 

None 

512 

0.098 

0.151 

Dynamics 

512 

0.098 

0.151 

Observations 

614 

0.117 

0.140 

Gain 

485 

0.093 

0.151 

All 

474 

0.090 

0.150 
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The  condition  number  for  the  Kalman  filter  covariance 
matrix  for  satellite  6633  is  lE+08  which  is  about  100  times 
larger  than  the  condition  number  for  satellite  10299.  This 
indicates  that  there  may  be  an  unobservable  state  for  this 
case  also.  An  attempt  was  made  to  process  the  data  using 
the  reduced  order  filter  described  previously  (where  all 
elements  but  ballistic  coefficient  are  estimated);  but  the 
filter  eventually  diverged  and  did  not  yield  meaningful 
results.  Prior  to  the  filter  divergence  the  condition 
number  was  the  same  as  before,  lE+08 ,  indicating  that  the 
unobservable  state (s)  is  not  ballistic  coefficient.  This 
is  reasonable  since  for  a  circular  orbit  the  argument  of 
perigee  and  mean  anomaly  are  undefined.  The  condition 
number  does  not  identify  the  unobservable  state(s);  however, 
the  physical  situation  indicates  that  the  argument  of  peri¬ 
gee  and  mean  anomaly  are  unobservable  for  satellite  6633. 

The  large  condition  number  can  be  used. to  explain  the 
large  residuals  for  the  recursive  filters.  A  large  condi¬ 
tion  number  indicates  numerical  instability  and  the  recursive 
filters  are  inherently  less  robust  than  the  least  squares 
filter.  This  is  further  supported  by  the  increase  in  the 
residuals  when  only  the  observation  bias  term  is  included 
in  the  Gauss  filter.  Normally,  the  inclusion  of  second 
order  terms  in  the  filter  should  improve  the  filter  per¬ 
formance,  or  at  worst,  not  change  the  filter  performance. 


In  this  case,  including  the  observation  bias  term,  the 
additional  term  in  the  filter  equations  caused  the  filter 
to  diverge,  resulting  in  larger  residuals. 

5.3.  Prediction  Errors 

The  filter  prediction  errors  are  those  errors  which 
result  from  propagating  the  filter  states  forward  in  time 
without  updating  the  state  vector  estimate  using  the  obser¬ 
vations.  The  prediction  errors  are  based  on  the  state 
estimate  at  time  t(i)  based  on  observations  up  to  time 
t(j),  where  t(j)  is  less  than  t(i).  The  state  vector  esti¬ 
mates  used  in  this  section  to  determine  the  prediction 
errors  are  the  final  values  of  the  state  vectors  used  in 
Section  5,2.  The  state  vector  estimate  at  time  t(j)  is  the 
state  vector  at  the  time  of  the  last  observation  processed 
in  the  estimation  problem.  This  estimate  is  then  propa¬ 
gated  forward  using  the  full  equations  of  motion  described 
in  Chapter  2.  The  prediction  errors  are  calculated  by 
comparing  the  predicted  value  of  the  observations  with 
actual  radar  observations  of  the  various  satellites. 

r  =  z  -  z  5.5 

—  — o  — c 

The  vector  of  residuals  is  made  up  of  residuals  in  slant 
range,  azimuth,  and  elevation. 

Figures  5-1  to  5-3  show  the  prediction  errors  for 
satellite  4507  (h=1000km,  e=0.015,  i=90  deg).  The  epoch 


values  of  the  state  vector  estimate  are  the  final  values 
of  the  state  vector  after  processing  three  days  of  obser¬ 
vational  data  using  the  full  order  filters;  the  filters 
estimate  the  classical  orbital  elements  plus  ballistic 
coefficient.  The  prediction  errors  are  evaluated  for  three 
days  following  the  last  measurement  update.  The  points 
shown  in  the  figures  are  the  average  residual  for  the 
previous  ten  observations. 

Figures  5-4  to  5-6  are  again  prediction  errors  for 
satellite  4507.  These  plots  are  based  on  the  state  vector 
estimate  from  the  reduced  order  filters;  ballistic  coeffi¬ 
cient  is  not  estimated.  For  both  sets  of  filters  the 
prediction  errors  for  the  two  recursive  filters,  the 
Extended  Kalman  filter  and  the  Gauss  Filter,  are  nearly 
identical  to  satellite  4507.  Also,  there  is  little  differ¬ 
ence  in  the  prediction  errors  between  the  original  filters 
and  the  reduced  order  filters.  Finally,  the  prediction 
errors  for  the  recursive  filters  appear  to  be  steady  while 
the  least  squares  filter  errors  appear  to  be  growing  line¬ 
arly.  This  linear  growth  in  the  prediction  errors  is  due 
to  an  error  in  the  estimate  of  the  semi-major  axis  at  the 
beginning  of  the  prediction  span.  This  initial  error  in 
the  semi-major  axis  results  in  a  timing  error  in  the  pre¬ 
diction  of  the  position  of  the  satellite.  In  other  words, 
the  predicted  position  of  the  satellite  may  be  very  accurate 


but  the  satellite  arrives  either  earlier  or  later  than 
the  predicted  time.  At  the  time  of  the  observation  for 
which  the  residuals  are  calculated  the  predicted  position 
of  the  satellite  is  calculated  for  a  different  point  in  the 
orbit.  The  position  of  the  satellite  in  the  orbit,  the 
true  anomaly  of  the  satellite,  is  calculated  by  integrating 
the  mean  motion  of  the  satellite  which  is  a  function  of 
the  semi -major  axis.  Therefore,  errors  in  the  semi-major 
axis  are  integrated  forward  with  time,  resulting  in  a  linear 
growth  in  the  predicted  true  anomaly  of  the  satellite. 

This  linear  growth  in  the  residuals  is  apparent  in  the 
plots  of  residuals  for  the  least  squares  filter. 

An  important  observation  can  be  made  from  comparing 
the  estimation  errors  with  the  prediction  errors  for  satel¬ 
lite  4507.  Even  though  the  estimation  errors  for  the  least 
squares  filter  are  smaller  than  the  estimation  errors  for 
both  the  Kalman  and  Gauss  filters,  the  prediction  errors 
for  the  Kalman  and  Gauss  filters  start  out  nearly  the 
same  as  the  least  squares  filter  and  remain  small  while  the 
prediction  errors  for  the  least  square  filters  have  a 
linear  growth.  This  indicates  that  even  though  the  least 
squares  filter  is  fitting  the  data  better  (smaller  estima¬ 
tion  residuals) ,  the  Kalman  and  Gauss  filters  have  a  better 
estimate  of  the  true  state  vector  (smaller  prediction 
residuals) .  The  estimation  residuals  are  obtained  by 


only  fitting  the  position  data  (range,  azimuth,  and  eleva¬ 
tion)  from  the  radar  observations.  The  propagation  of  the 
orbit  forward  in  time  without  measurement  update  requires 
a  knowledge  of  the  velocity  as  well  as  the  position  of 
the  satellite.  Therefore,  the  least  squares  filter 
appears  to  be  providing  a  better  estimate  of  the  position 
of  the  satellite,  and  the^  Kalman  and  Gauss  filters  appear 
to  be  providing  poorer  estimates  of  the  position  but  better 
estimates  of  the  velocity  of  the  satellite. 

Similar  results  were  obtained  for  the  prediction  errors 
for  satellite  4507  when  the  state  estimates  from  the  filters 
at  different  times  were  used  to  generate  the  prediction 
residuals . 

The  prediction  residuals  for  satellite  10299  (h  =  250kim, 
e=0.01,  i=72  deg)  are  shown  in  Figures  5-7  to  5-9.  The 
filter  state  vector  estimates  are  propagated  forward  one 
day.  The  points  on  the  plots  represent  the  average  of  the 
five  previous  residuals.  The  state  estimate  for  the  Gauss 
filter  is  the  state  estimate  using  all  three  biases; 
dynaunics,  observations,  and  gain,  in  the  filter  equations. 

The  linear  growth  in  the  residuals  is  present  for  the  least 
squares  filter  and  the  Kalman  Filter  for  this  satellite; 
however,  the  residuals  for  the  Gauss  filter  are  nearly 
constant.  Since  the  orbit  of  the  satellite  is  affected  by 
drag  the  growth  in  the  residuals  due  to  a  timing  error  also 
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has  a  quadratic  term.  In  addition  to  the  integration  of  the 
initial  errors  in  the  semi-major  axis,  errors  in  the  semi¬ 
major  axis  and  ballistic  coefficient  result  in  errors  in  the 
magnitude  of  the  drag  perturbation.  These  errors  are  inte¬ 
grated  once  to  obtain  the  semi-major  axis  and  once  again  to 
obtain  the  true  anomaly,  resulting  in  a  quadratic  growth  in 
the  timing  error  and  therefore  in  the  prediction  residuals. 

Figures  5-7  through  5-9  clearly  show  the  advantage  of 
including  the  second  order  terms  in  the  filter  equations. 

The  predicted  orbit  based  on  the  state  vector  estimate 
from  the  Gauss  filter  provides  a  much  better  match  to  the 
observed  orbit  (position  and  velocity)  than  the  estimate 
from  either  the  least  squares  filter  or  the  Kalman  filter. 

Figures  5-10  through  5-12  are  the  prediction  residuals 
for  satellite  6633  (h=300km,  e=0.0003,  i=50  deg).  The  state 
estimates  are  propagated  forward  for  two  days  and  the  points 
on  the  curves  represent  the  average  of  the  previous  five 
residuals.  Again,  the  state  vector  estimate  for  the  Gauss 
filter  is  based  on  the  estimate  using  all  three  bias  terms 
in  the  filter  equations.  Significant  growth  in  the  resid¬ 
uals  appears  in  the  plots  for  the  Kalman  and  Gauss  filters; 
however,  the  prediction  errors  for  the  least  squares  filter 
are  nearly  constant.  The  initial  errors  for  the  Kalman 
and  Gauss  filters  are  also  much  larger  than  for  the 
previous  cases.  The  large  initial  residuals  in  the 
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propagated  states  are  indicative  of  large  errors  in  the 
state  vector  estimate  at  the  end  of  the  estimation 
sequence.  This  result  is  not  unexpected  since  the  condi¬ 
tion  number  for  the  covariance  matrix  for  satellite  6633 
is  100  times  larger  than  that  for  satellite  10299.  This 
large  condition  number  indicates  one  or  more  unobservable 
states.  Loss  of  numerical  accuracy  and  instability  can  be 
expected  when  the  condition  number  is  large.  Physically 
this  is  reasonable  because  the  orbit  for  satellite  6633  is 
nearly  circular  making  it  very  difficult  to  determine  argu 
ment  of  perigee  and  mean  anomaly.  In  spite  of  these 
numerical  problems,  the  estimate  from  the  least  squares 
filter  does  not  diverge  and  provides  an  adequate  estimate 
of  the  orbit  of  the  satellite  for  re-acquisition  based  on 
the  predicted  look  angles. 

5.4.  Computer  Time  of  Operation 

The  computer  time  of  operation  is  the  computer  execu¬ 
tion  time  required  to  process  the  data  used  in  Section  5.2 
to  determine  the  estimation  residuals.  The  data  are 
described  in  Chapter  4.  The  time  of  operation  for  the 
Gauss  filter  is  presented  for  several  different  cases. 
GSF(D)  is  the  Gauss  filter  with  the  dynamic  bias  only, 

GSF{0)  if  the  Gauss  filter  with  only  the  observation  bias, 

GSF(G)  is  the  Gauss  filter  with  the  gain  bias  only,  and 
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GSF (A)  is  the  Gauss  filter  with  all  of  the  bias  terms 
included . 

The  computer  times  for  the  filters  for  each  satellite 
case  are  shown  in  Table  5-8- 

Table  5-8 


Computer  Time  of 

Operation 

(sec ) 

Filter 

Satellite 

1 

Test  Case 

4507 

10299 

6633 

LSF 

89.3 

125.3 

98.2 

216.1 

EFK 

18.1 

39.0 

19.5 

29.6 

GSF(D) 

- 

- 

’32.9 

55.7 

GSF(O) 

- 

- 

25.8 

36.4 

GSF(G) 

- 

- 

26.0 

36.7 

GSF (A) 

32.1 

77.4 

38.8 

59.6 

The  Air  Force  Space  Command  is  responsible  for  main¬ 
taining  current  orbital  element  sets  for  about  6000 
objects  orbiting  the  Earth.  Of  the  objects,  300  -  350 
are  active  satellites  and  the  remainder  are  debris.  The 
reduction  in  the  amount  of  computer  time  required  for  the 
problem  of  element  set  maintenance  for  the  satellites 
listed  in  Table  5-8  becomes  even  more  significant  when 
applied  to  a  large  number  of  satellites. 


5.5.  Computer  Storage  Requirements 

The  computer  storage  requirements  for  the  three  filter 
are  compared  in  this  section.  The  data  shown  in  Table  5-9 
are  the  computer  memory  required  to  load  and  run  the  object 
codes  and  supporting  software  for  each  of  the  filters. 

There  is  not  a  significant  difference;  however,  the  data 
do  not  include  the  memory  required  for  storing  data  prior 
to  batch  processing  for  the  least  squares  filter.  This 
amount  of  storage  can  become  significant  for  element  set 
maintenance  for  a  large  number  of  satellites. 

Table  5-9 

Computer  Storage  Requirements  (k  of  memory) 


Filter 


Storage  (k) 


CHAPTER  6 


CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 

6.1.  Cone lusions 

Several  conclusions  can  be  drawn  from  the  results  pre¬ 
sented  in  Chapter  5. 

1.  The  state  vector  estimate  from  the  Gauss  filter 
provides  the  best  estimate  for  propagating  the  state  vector 
forward  in  time  for  satellites  4507  and  10299.  The  predic¬ 
tion  residuals  for  satellite  10299  clearly  show  that  the 
smallest  residuals  result  from  the  propagation  of  the  GSF 
state  vector.  The  orbit  of  satellite  10299  (h=250]an)  is 
typical  of  a  large  number  of  satellites  currently  in  orbit 
around  the  Earth.  These  orbits  are  non-circular  and  have 
low  altitudes  where  the  perturbations  due  to  the  non- 
spherical  Earth  and  to  atmospheric  drag  are  most  significant. 
For  satellite  4507  the  prediction  residuals  for  the  EKF  and 
GSF  are  nearly  identical  and  smaller  than  the  LSF  residuals. 
The  nearly  identical  results  from  the  Kalman  and  Gauss 
filters  for  satellite  4507  are  not  unexpected.  The  orbit 
of  satellite  4507  is  at  a  high  altitude  (h=1000km)  where 
the  perturbations  due  to  a  non-spherical  Earth  and  atmos¬ 
pheric  drag  are  nearly  zero;  therefore,  the  second  order 
terms  in  the  Taylor  series  should  have  little  effect  on 
the  problem  of  orbit  determination.  The  estimates  from 


the  Kalman  and  Gauss  filters  for  satellite  6633  diverged 
due  to  singularities  in  the  orbit  determination  problem 
for  this  orbit  (e=0.0003);  however,  the  least  squares 
filter  provided  an  adequate  estimate  of  the  state  vector 
for  predicting  the  future  position  and  velocity  of  the 
satellite. 

2.  The  least  squares  filter  estimate  of  the  state 
vector  results  in  the  smallest  estimation  residuals  when 
processing  data  for  all  three  satellites.  The  largest 
difference  in  estimation  residuals  is  for  satellite  6633 
which  is  in  a  nearly  circular  orbit.  The  problem  of 
determining  the  classical  orbital  elements  of  a  circular 
satellite  is  singular.  For  nearly  circular  orbits,  the 
Kalman  and  Gauss  filters  experience  numerical  problems  and 
the  estimates  diverge.  The  LSF  appears  to  be  more  robust 
and  less  susceptible  to  numerical  problems. 

3.  The  Gauss. filter  estimation  residuals  are  equal  to 
or  smaller  than  the  estimation  residuals  for  the  Kalman 
Filter.  The  biggest  improvement  is  for  satellite  10299. 

The  orbit  of  satellite  10299  is  low  altitude  (h=250  km) 

and  noncircular  (e=0.01)  representing  a  case  where  the 
determination  of  the  classical  orbital  elements  is  non¬ 
singular.  In  addition,  the  perturbing  accelerations  due  to 
atmospheric  drag  and  the  non-spherical  Earth  are  significant. 
Little  or  no  improvement  in  the  estimation  residuals  for 
the  Gauss  filter  is  observed  for  satellite  4507  where  the 
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perturbing  accelerations  are  nearly  zero.  This  is  in  agree¬ 
ment  with  the  conclusions  of  Tapley  and  Choe  (5)  and 
Athens  et  al .  (7).  They  both  report  that  the  second  order 
filters  provide  the  greatest  improvement  over  the  Kalman 
filter  when  the  perturbations  are  the  most  significant. 
Little  can  be  concluded  from  the  results  for  satellite 
6633  due  to  the  singularities  in  the  problem. 

4.  The  Extended  Kalman  filter  is  much  faster  than  the 
least  squares  filter  and  the  Gauss  filter.  The  EKF  is 
about  twice  as  fast  as  the  GSF  in  processing  the  data  for 
all  three  satellites  and  nearly  five  times  as  fast  as  the 
LSF.  The  GSF  processes  the  data  two  to  three  times  faster 
than  the  LSF.  These  differences  in  process  time  become 
even  more  significant  when  applied  to  the  problem  of 
maintaining  current  orbital  element  sets  for  6000  satel¬ 
lites. 

5.  The  computer  storage  requirements  are  nearly  the 
same  for  all  three  filters.  The  computer  storage  require¬ 
ments  do  not  address  the  need  to  store  the  observations  of 
the  satellites  prior  to  batch  processing  by  the  least 
squares  filter. 

6.  The  most  significant  second  order  terms  in  the 
filter  equations  are  found  in  the  observation  equations. 

The  inclusion  of  the  bias  terms  for  the  dynamics  and  gain 
in  the  filter  equations  result  in  no  measurable  improvement 


in  the  Gauss  filter  performance.  Taylor  (6)  also  reported 
that  greatest  effect  due  to  the  nonlinearities  in  the 
orbit  determination  problem  are  caused  by  nonlinearities 
in  the  observation  equations.  This  is  in  contrast  to  the 
work  of  Tapley  and  Choe  (5),  who  reported  that  the  most 
significant  bias  term  in  the  interplanetary  problem  is 
the  gain  bias  term  while  Athens,  et  al.  (7)  reported  that 
the  dynamic  bias  term  was  the  most  significant  bias  term 
for  a  vertically  falling  body.  This  problem  can  be  ex¬ 
plained  by  the  fact  that  these  two  estimation  problems 
are  qualitatively  quite  different  from  the  determination 
of  the  orbit  of  an  Earth  satellite. 

The  inclusion  of  the  bias  term  for  the  observations 
in  the  filter  improves  the  Gauss  filter  performance  for 
satellite  10299,  but  results  in  no  change  for  satellite 
4507.  Little  can  be  concluded  for  satellite  6633  due  to 
the  singularity  in  the  problem  of  determining  the  classi¬ 
cal  orbital  elements  for  a  nearly  circular  orbit. 

7.  The  reduced  order  filters  (in  which  ballistic 
coefficient  was  removed  from  the  estimation  state  vector) 
for  the  high  altitude  satellite  4507  reduce  the  estimation 
residuals  for  the  Kalman  and  Gauss  filters. 

6.2.  Recommendations  for  Future  Work 

1.  The  orbit  determination  problem  for  satellite 
6633  presents  several  problems  due  to  the  singularities 
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in  the  classical  orbital  elements.  Also,  estimating  the 
osculating  elements  results  in  large  short  period  varia¬ 
tions  in  the  state  vector  over  relatively  short  time 
periods.  This  presents  two  problems.  First,  the  step 
size  for  the  numerical  integrator  is  a  function  of  the 
magnitude  of  the  variations  of  the  state  vector.  Also,  the 
assumptions  made  in  equation  3.13  to  calculate  the  state 
transition  matrix  introduces  errors  into  the  problem  when 
the  elements  of  the  state  vector  are  not  constant  over  the 
period  of  interest.  A  transformation  of  variables  from  the 
state  vector  used  in  this  research  to  a  set  of  mean, 
non-singular  elements  should  improve  the  filter  perform¬ 
ance.  One  set  of  elements  well  suited  to  the  orbit 
determination  problem  is  the  set  of  mean  equinoctial 
elements  (3) .  The  equations  of  motion  and  observation 
equations  described  in  Chapter  2  should  be  converted  to  the 
mean,  equinoctial  set  of  orbital  elements  and  the  filter 
comparisons  described  in  Chapter  5  should  be  repeated. 

2.  A  large  portion  of  the  computer  operation  time 
is  required  to  propagate  the  state  vector  from  one  obser¬ 
vation  to  the  next.  This  method  of  special  perturbations 
is  accurate  but  time  consuming.  Semi-analytical  orbit 
theory  reduces  the  amount  of  computer  time  required  for 
numerical  integration  of  the  equations  of  motion  by 
solving  a  portion  of  the  equations  of  motion  analytically 


(6  and  21) .  The  solution  of  the  filter  equations  of 
motion  should  include  semi-analytic  orbit  theory  where 
applicable. 

3.  Third-body  perturbations  due  to  the  sun  and  moon 
should  be  included  in  the  equations  of  motion  to  solve  the 
orbit  determination  problem  for  high  altitude  satellites. 

4.  Different  methods  of  calculating  the  state  transi¬ 
tion  matrix  should  be  investigated.  This  should  be  done 

in  conjunction  with  suggestion  one.  Estimating  mean  orbi¬ 
tal  elements  instead  of  osculating  orbital  elements  may 
make  the  assumption  in  equation  3.13  more  valid. 

5.  The  assumption  that  the  process  noise  matrix, 

Q,  and  the  measurement  noise  matrix,  R,  are  diagonal 
matrices  should  be  relaxed.  The  effect  of  off  diagonal 
elements  in  both  matrices  should  be  investigated. 

6.  Additional  tracking  data  for  low,  non-circular 
orbits  should  be  processed  to  confirm  and  quantify  the 
precision  and  advantages  of  the  Gauss  filter. 
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APPENDIX  A 


FUNCTIONS  AND  COEFFICIENTS  FOR  TESSERAL  HARMONICS 


This  appendix  contains  the  expressions  for  the  incli¬ 
nation  functions  and  the  eccentricity  functions  found  in 
the  tesseral  terms  in  the  geopotential  as  a  function  of  the 

I 

classical  orbital  elements,  equation  2.17.  Also  contained 
in  this  appendix  are  the  values  for  the  qeopotential 


coefficients  S.  and  C.  , 
2,m  x,m 


found  in  equation  2.17. 


Table  A-1 

Inclination  Functions  F. _  (13) 

Jimp 

220  3(l+cosi)^/4 

221  3sin^i/2 

2  2  2  3  (l-cosi)V4 

310  -ISsin^i (1+cosi) /16 

311  ISsin^i (l+3cosi) /16-3 (1+cosi) /4 

312  ISsin^i (l-3cosi) /16-3 (1-cosi) /4 

313  -ISsin^i (1-cosi) /16 

320  15sini(l+cosi) ^/8 

321  ISsini (l-2cosi-3cos^i) /8 

322  -15sini(l-2cosi-3cos^i)/8 
-ISsini (1-cosi) ^/8 


3 


2 


3 


83 


Table  A-1  (Continued) 

331  4 Ssin^i (1+cosi) /8 

3  3  2  4  Ssin^i  (1-cosi) /8 

3  3  3  ISd-cosi)  ^/8 

410  -35sin^ i (1+cosi) /32 

^  4  11  35sin^i  (l+2cosi) /16-15  (1+cosi)  sini/8 

412  cosi (15sini/4-105sin^ i/16) 

413  -35sin^i (l-2cosi) /16+15sini (l-cosi) /8 

414  35sin^i (1-cosi) /32 

4  2  0  -lOSsin^i  (1+cosi)  V32 

421  105 (sin^icosi (1+cosi) /8-15 (1+cosi) ^/8 

422  105sin^i(l-3cos^i)/16+15sin^i/4 

4  2  3  -lOSsin^icosi  (l-cosi]/8-15  (l-c6si)  ^/8 

4  24  -105sin^i(l-cosi) ^/32 

430  lOSsini (l+cosi) ^/16  • 

431  105sini(l-3cos^i-2cos^i)/8 

432  -315sin^icosi/8 

433  -lOSsini (l-3cos^i+2cos^i) /8 

434  -lOSsini (1-cosi) ^/1 6 

4  4  0  105(l+cosi)  "^/le 

441  lOSsin^i (1+cosi) ^/4 

442  315sin^i/8 

443  lOSsin^i ( 1-cosi) ^/4 

444  105 (1-cosi) ^/16 
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Table  A- 2 

Eccentricity  Functions  (13) 

1  p  q  1  P  q  G,pg(e) 

2  0  -1  2  2  1  -e/2 

2  0  0  2  2  0  1 

2  0  1  2  2  -1  7e/2 

2  1-12  11  3e/2 

2  10  1 

3  0  -1  3  3  1  -e 

3  0  0  3  3  0  1 

3  0  1  3  3  -1  5e 

3  1  -1  3  2  1  e 

3  1  0  3  2  0  1 

3  1  1  3  2  -1  3e 

4  0  -1.4  4  1  •  -3e/2 

4  ,  0  0  4  4  0  1 

4  0  1  4  4  -1  13e/2 

4  1-14  3  1  e/2 

4  1  0  4  3  0  1 

4  114  3-1  9e/2 

4  2  -1  4  2  1  5e/2 


4 


2 


0 


1 


85 


Table  A-3 


Tesseral 

Coefficients 

(22) 

1 

m 

S.  E-06 
Im 

2 

2 

1.5577 

-.8806 

3 

1 

2.1277 

.4157 

3 

2 

.3047 

-12168 

3 

3 

.0957 

.1995 

4 

1 

-.5027 

-.4627 

4 

2 

.0738 

.1579 

4 

3 

-.0591 

-.0092 

4 

4 

-.0017 

.0072 
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APPENDIX  B 

STATE  TRANSITION  MATRIX 

The  state  transition  matrix,  ,  is  used  to  calculate 
the  linearized  observation  matrix,  H.  Specifically,  it  is 
the  matrix  of  partial  derivatives  of  the  state  vector  at 
some  time  t  with  respect  to  the  state  vector  at  epoch 
time,  t^.  The  method  of  determining  the  state  transition 
matrix  in  this  research  is  based  on  eq.  3.13: 


$  »  I  +  (t-t^)?,  B-1 

where: 

I  is  the  identity  matrix,  and 
F  is  the  linearized  dynamics  matrix. 

The  linearized  dynamics  matrix  used  in  the  filter  equa¬ 
tions  in  Chapter  3  and  in  equation  B-1  to  calculate  the  state 
transition  matrix  includes  the  two-body  acceleration,  varia¬ 
tions  due  to  atmospheric  drag,  and  secular  variations  due 
to  the  zonal  terms  in  the  geopotential.  The  terms  in  the 
F  matrix  due  to  two-body  acceleration,  the  tangential  com¬ 
ponent  of  atmospheric  drag,  and  secular  variations  due  to 
J2  are  presented  in  the  following  equations.  All  elements 
of  the  F  matrix  not  shown  below  are  zero  using  this  dynamics 
model . 


F(l,l)  =  -nQBaOpexp{-z)  [z(Ij^-Iq)+Iq/2]  B-2 

F(l,2)  »  nQBza^Ppexp(-z)  [lQ-Ij^]/e  B-3 
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F(l,7)  »  -nQa^Ppexp (-2) B_4 

F(2,l)  »  -nQBppexp(-z)  l2(I^-I^) B-5 

F(2,2)  =  nQBappexp(-z) [z +1^] /e  B-6 

F(2,7)  =  -nQappexp  (-2)  B-7' 

F(4,l)  =  (21/4) J2(n/a(R/p) ^cos(i)  B-8 

F(4,2)  *  -6nJ2 (R/p) ^ecos (i) / (i-e^)  B-9 

F(4,3)  =  1. 5nJ2(R/p) ^sin(i)  B-10 

F(5,l)  =  (21/4) J2(n/a) (R/p)^{2.5sin^(i)-2)  B-11 

F(5,2)  =  -6nJ2(R/p)^(2.5sin^(i)-2)/(l-e^)  B-12 

F(5,3)  =  -7.5nJ2(R/p) ^sin(i)cos (i)  B-13 

F(6,l)  =  -1.5n/a+(2l/4)J2(n/a) (R/p)2 

(1. 5sin2(i) -1) (i-e^)^/^  B-14 


F(6,2)  =  4.5nJ2(R/p)^(1.5sin^(i)-l) (e/(l-e^)^/^)  B-15 
F(6,3)  =  -4 . 5nJ2  (R/p)  ^sin(  i)  cos  ( i)  (1-e^) 


B-16 


APPENDIX  C 


SECOND  ORDER  TERMS  FOR  FILTER  EQUATIONS 

The  second  order  filters  described  in  Section  3-4 
include  terms  which  require  the  calculation  of  the  second 
partial  derivatives  of  the  dynamics  equations  with  respect 
to  the  state  vector  at  epoch  and  the  second  partial  deriva 
tives  of  the  observations  with  respect  to  the  state  vector 
at  epoch.  This  appendix  describes  the  approach  used  to 
evaluate  these  derivatives. 

The  second  partial  derivatives  of  the  dynamics 
equations,  f(x,t),  with  respect  to  the  state  vector  at 
epoch,  x(t^) ,  used  to  evaluate  the  dynamic  bias  term, 
will  be  described  first.  The  final  result  is  a  set  of 
seven,  7x7  matrices.  These  seven  matrices  are  evaluated 
one  at  a  time  using  the  rows  from  the  linearized  dynamics 
matrix,  F,  described  in  Appendix  B.  The  only  accelera¬ 
tions  needed  to  calculate  the  second  partial  derivatives 
are  the  two-body  acceleration,  the  tangential  component 
of  the  drag  acceleration,  and  the  secular  component  of 
the  acceleration  due  to  oblateness,  (J2) .  The  first  of 
the  seven  matrices,  the  second  partial  derivative  of  the 
variation  in  the  semi-major  axis  with  respect  to  the 
state  vector  at  epoch,  is  evaluated  using  row  one  of  the 
F-matrix.  Similarly,  the  second  partial  derivative  of  the 
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variation  in  eccentricity  with  respect  to  the  epoch  state 
vector  is  evaluated  using  the  second  row  of  the  F-matrix. 

It  can  be  immediately  observed  that  the  matrices 
corresponding  to  the  variations  in  inclination  and  ballis¬ 
tic  coefficient  are  the  null  matrices;  therefore,  the  third 
and  seventh  elements  of  the  dynamic  bias  vector  will  be 
zero. 

The  non-zero  elements  of  the  seven  matrices  which 
represent  the  second  partial  derivatives  of  the  dynamics 
with  respect  to  the  state  vector  at  each  epoch  time  are 
listed  below. 

Matrix  due  to  variation  in  semi-major  axis: 

Pl(l,l)  *  nQBp  exp(-z)  I2.0z^(I. -I  )+I  (z+1/4)] ,  C-1 

p  loo 


Fid, 2)  =  -nQBa/ep  exp(-z)  [z(I, -I  ) +I  e/21,  C-2 

p  loo 

Fl(l,7)  =  -nQaPp  exp(-z)  [z(Ij_-I^)+I^/2] ,  C-3 

Fl(2,l)  *  Fid, 2),  C-4 

Fl(2,2)  »  nQBz(a/e)2p  exp(-z)I,,  C-5 

Fl(2,7)  =  nQza^Ppexp(-z)  (I^-Ij^)/e,  C-6 

Fl(7,l)  =  Fid, 7)  ,  C-7 

Fl(7,2)  =  Fl(2,7) .  C-8 
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Matrix  due  to  variation  in  eccentricity; 

F2(l,l)  =  (n/a)QBppexp(-z) [z(22+3) C-9 


F2(l,2)  =  nQBppexp(-z)  z(z-l)  (I^-Ij^l/e,  C-10 

F2(l,7)  *  nQPpexp(-z)  [z(I^-I^)+1.5I^] ,  C-11 

F2(2,l)  =  F2(l,2)  ,  C-12 

ft 

F2(2,2)  =  nQBappexp(-z)  r2z^(I^-Ij^) 

+z{I^-2lj^)-IJ/e^,  C-13 

F2(2,7)  »  nQ(a/e)Ppexp(-z)  t2(I^-Ij^)+I^]  ,  C-14 

F2(7,l)  »  F2(l,7) ,  C-15 

F2(7,2)  «  F2(2,7) .  C-16 

Matrix  due  to  variation  in  inclination: 

F3 (I,J)  =  0.0.  C-17 

Matrix  due  to  variation  in  ascending  node: 

F4(l,l)  »  -(189/9) J2(n/a^) (R/p) ^cos(i) ,  C-18 

F4{1,2)  -  21J2(n/a) (R/p)^ecos(i)/(i-e^) ,  C-19 

,  F4  (1,3)  =  -(21/4) J2(n/a) (R/p) ^sin(i) ,  C-20 

F2(2,l)  »  F4 (1,2) ,  C-21 
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F4(2,2) 

=  -6J2n(R/p) ^cos(i) (l+5e^) / (1-e^) 

C-22 

F4(2,3) 

*  6nJ2 (R/p) ^esin (i) / (1-e^) , 

C-23 

F4(3,l) 

=  F4(l,3) , 

C-24 

F4(3,2) 

*  F4(2,3) , 

C-25 

F4(3,3) 

=  1.5nJ2(R/p) ^cos(i) . 

C-26 

Matrix  due  to  variation  in  argument  of  perigee: 

F5(l,l)  *  -(189/8)J2(n/a^) (R/p)^[2.5sin2(i)-2] , 

C-27 

F5(l,2) 

=  21(n/a)  J2(R/p)^e[2.5sin2(i)-2]/(l-e2)  , 

C-28 

F5(l,3) 

®  (105/4)  J2(n/a).(R/p)  ^sin(i)cos  (i)  , 

C-29 

F5(2,l) 

»  F5(l,2)  , 

C-30 

F5(2,2) 

=  -6nJ2(R/p)^[2..5sin^(i)-2]  (l+5e^)  /  (1-e^) 

C-31 

F5(2,3) 

=  -30nJ2 (R/p) ^esin (i) cos (i) / (l-e^) , 

C-32 

F5(3,l) 

=  F5(l,3) , 

C-33 

F5(3 ,2) 

=  F5(2,3) , 

C-34 

F5(3,3) 

-  -7.5nJ2(R/p)^(cos^(i)-sin^(i)]  . 

C-3  5 

Matrix  due  to  variation  in  mean  anomaly: 

F6(l,l)  -  (15/4)n/a-{189/8) JjCn/a^) (R/p)^ 

[1.5sin^(i)-llx(i-e^)^/^,  C-36 


C-37 


F6(l,2) 

F6(l,3) 

F6(2,l) 

F6(2,2) 

F6(2,3) 
F6(3,l) 
F6(3,2) 
F6(3  ,3) 


-(63/4)  (n/a)J2(R/p)  ^(1.5sin^(i)-l]e/ 

( 63/4 ) (n/a) J2 (R/p) ^sin ( i) cos ( i) 

F6(l,2),  C-39 


-4.5nJ2(R/p) ^(1.5sin^{i)-l) {l+4e^)/ 

(l-e2)3/2^  C-40 

-13.5nJ2(R/p)  ^sin(i)cos  (i)  e/d-e^)^'^^,  C-41 

F6(l,3),  C-42 

F6(2,3),  C-43 

4.5nJ2(R/p)^(l-e2)^/2jg^j^2(^j^^g2^^jj  ^ 


Matrix  due  to  variation  in  ballistic  coefficient: 


F7(I,J)  =  0.0.  C-45 

The  other  second  order  terms  in  the  Gauss  filter 
described  in  Section  3.4  require  the  calculation  of  the 
second  partial  derivatives  of  the  observations  with  respect 
to  the  epoch  state  vector.  This  results  in  three  7x7 
matrices,  corresponding  to  slant  range,  azimuth,  and 
elevation.  These  three  matrices  are  used  to  calculate 
the  bias  terms  in  the  observations  and  the  filter  gain. 
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The  starting  point  for  calculating  the  three  matrices  is 
the  linearized  observation  matrix,  H.  The  dynamics  model 
described  above  is  used  once  again.  The  evaluation  of 
the  second  partial  derivatives  would  be  relatively  straight 
forward;  however,  the  elements  of  the  H-matrix  are  not 
available.  Instead,  the  H-matrix  is  expressed  as  a  pro¬ 
duct  of  five  matrices: 

H  =  Ml  M2  M3  M4  M5.  C-46 

The  approach  used  is  to  first  evaluate  the  partial 
derivatives  of  H  with  respect  to  each  of  the  epoch  states. 
This  results  in  seven  3x7  matrices.  The  partial  derivatives 
are  evaluated  using  the  product  rule  of  differentiation 
and  the  chain  rule.  The  first  two  matrices  in  equation 
C-46  are  not  functions  of  the  epoch  state  vector;  however, 
the  last  three  matrices  are  functions  of  the  epoch  state 
vector.  The  derivatives  of  the  last  three  matrices  in 
equation  C-46  with  respect  to  the  elements  of  the  state 
vector  are  presented  below.  The  notation  used  in  the  fol¬ 
lowing  equations  is:  M3A  is  the  partial  derivative  of 
matrix  M3  with  respect  to  the  semi-major  axis,  etc. 


M3A(I,J)  =*  0.0,  1=1,... 3;  J=l,..., 


7 


C-47 
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M3 1 


'sinfi sinus ini  -cosfi sinus ini 
sinflcosusini  -cosficosusini 
sinucosi  -cosucosi 


sinucosi 

cosucosi 

-sini 


C-49 


M3n 


-sinOsosu-cosA  .sinucosi  cosficosu-sirfi  sinucosi  0 

sirflsinu-co^cosucosi  -cosnsinu-sinfi  sinucosi  0 

cosfisini  sirfisini  0 


C-50 


M3  CO 


-cosfisinu-sinficoaicosi  -sinfisinu+(X)Sficosucosi  oosusini 
-cosfisinu+sirfi sinucosi  -sirficosu-co^ sinucosi  -sinusini 


C-51 


0 


0 


0 


M3M  =  (a/r)  ^(l-e^)*/^M3co.  C^52 
M4A(1,2)  =  (r/a-l)/e,  C-53 
M4A(1,6)  =  esin(f)/(l-e^)  C-54 
M4A(2,4)  =  (r/a)cos(i),  C-55 
M4A(2,5)  =  r/a  C-56 
M4A(2,6)  =  (a/r) (l-e^)l/2,  C-57 
M4A(3,3)  =  (r/a)sin(u),  C-58 
M4A{3,4)  =  - (r/a) sin(i)cos (u) .  C-59 
M4E(1,1)  =  k/a,  C-60 
M4E(1,2)  =  k/e-(r-a)/e^,  C-61 
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M4E(1,6)  =  asin(f)/(l-e^)^^^,  C-62 

M4E{2,4)  =  kcos(i),  C-63 

M4E(2,5)  =  k,  C-64 

M4E(2,6)  =  a  (cos  (f ) +ae/r)  /  (l-e^)  C-65 

M4E(3,3)  =  ksin(u),  C-66 

M4E(3,4)  =  -ksin  (i)  cos  (u)  ,  C-67 

where: 


k  =  3r/3e  = 

=  -a (2e+e^cos (f ) +cos (f ) ) / (1+ecos (f ) ) ^ . 

C-68 

M4I(2,4)  = 

-rsin (i) , 

C-69 

M4I(3,4)  = 

-rcos (i)cos (u) . 

C-7  0 

M4S^(I,J)  = 

0.0. 

C-71 

M4(o(3,3)  = 

rcos (u) , 

C-72 

M4a)(3,4)  = 

rsin ( i) sin (u) . 

C-73 

M4M(1,1)  = 

k/a. 

C-74 

M4M(1,2)  = 

k/e. 

C-75 

M4M(1,6)  = 

2 

(a/r)  aecos(f). 

C-76 

M4M(2, 4)  =  kcos (i) 


C-77 


M4M(2,5)  =  k. 


C-78 


M4M(2,6)  =  -aesin(f),  C-79 

M4M(3,3)  =  ksin  (u) -r  (a/r)  ^  (1-e^)  (u)  ,  C-80 

M4iM(3,4)  =  -ksin  ( i)  cos  (u) +r  (a/r)  ^ 

(1-e^) ^^^sin(i) sin{u) ,  C-81 

Where : 

k  =  3r/3M  =  resin (f) / (1+ecos (f) )  .  C-82 

M5A(I,J)  =  F1(I,J) (t-t^) .  C-83 

M5E(I,J)  =  F2(I,J) (t-t^) .  C-84 

M5I(I,J)  =  F3 (I,J) (t-t^)  =  0.  C-85 

M5fi(I,J)  =  F4(I,J) (t-t^) .  -  C-86 

M5w(I,J)  =  F5(I,J) (t-t^) .  C-87 

M6M(I,J)  =  F6(I,J) (t-t^) .  C-88 

M6B{I,J)  =  F7 (I,J) =  0  C-89 


The  above  equations  are  combined  to  form  the  seven  3x7 
matrices  representing  the  second  partial  derivatives  of  the 
observations  with  respect  to  each  state. 

3^h/3a^  =  Ml  M2[M3A  M4  M5+M3  M4A  M5+M3  M4  MSA], 


C-90 


97 


2  2 
3‘^h/3e 

s 

Ml 

M2 [M3E 

M4 

M5+M3 

M4E 

M5+M3 

M4 

M5E]  , 

C-91 

3^h/3i^ 

= 

Ml 

M2 [M3 1 

M4 

M5+M3 

M4I 

M5+M3 

M4 

MSI]  , 

C-92 

2  2 

3  h/3n 

= 

Ml 

M2  [M3n 

M4 

M5+M3 

M4n 

M5+M3 

M4 

MSn]  , 

C-93 

2  2 
3"^h/3<u 

= 

Ml 

M2  [M3a) 

M4 

M5+M3 

M4uj 

M5+M3 

M4 

MSoj]  , 

C-94 

2  2 
3^h/3M 

s 

Ml 

M2  [M3M 

M4 

M5+M3 

M4M 

M5+M3 

M4 

M5M]  , 

C-9S 

2  2 
3  h/3B‘' 

— 

Ml 

M2[M3B 

M4 

M5+M3 

M4B 

M5+M3 

M4 

MSB]  . 

C-96 

Once  the  seven  3x7  matrices  described  above  have  been 
calculated,  the  matrices  required  for  the  calculation  of 
the  biases  in  the  gain  and  observations  can  be  determined. 
The  matrix  corresponding  to  the  slant  range  is  determined 
by  combining  the  first  rows  of  the  seven  3x7  matrices  into 
a  7x7  matrix.  The  first  row  is  the  first  row  of  equation 
C-90,  the  second  row  is  the  first  row  of  equation  C-91, 
etc.  The  matrix  corresponding  to  the  azimuth  is  formulated 
by  using  the  second  rows  of  the  seven  matrices  described 
above.  Finally,  the  matrix  corresponding  to  elevation  is 
formulated  by  combining  the  last  rows  of  the  seven  matrices 


described  above. 
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APPENDIX  D 
RADAR  SITE  LOCATIONS 

The  locations  of  the  radar  sites  are  provided  by 
SPACECOM/DOA.  The  components  of  the  site  position  vectors 
are  defined  in  terms  of  a  coordinate  frame  fixed  to  the 
Earth  with  the  origin  at  the  center  of  the  Earth.  The 
principal  direction,  x,  is  the  Greenwich  meridian,  ^ 
is  aligned  with  the  North  Pole,  and  y  lies  in  the  equator¬ 
ial  plane  forming  a  right  handed  coordinate  frame.  The 
radar  site  locations  used  in  this  research  are  shown  in 
Table  D-1. 


Table  D-1 

Radar  Site  Locations 


Site  Number 

X  (kro) 

Y  (km) 

Z  (km) 

354 

6119.393 

-1517.496 

-871.566 

359 

-2382.980  - 

-1420.937 

5724.060 

363 

2881.604 

-5372.517 

1868.026 

393 

-3849.481 

398.421 

5052.965 

396 

-579.421 

-4175.712 

4770.682 

399 

362.839 

-5484.293 

3225.187 
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